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1.1.  Objective  of  the  Dissertation 

The  objective  of  this  dissertation  is  to  provide  tables  for 
time-dependent  expected  server  load  in  M/G/l  queueing  systems.^  The 
results  are  presented  in  a form  that  can  be  applied  by  practitioners 
involved  in  the  design  and  control  of  operating  systems.  This  presenta- 
tion attempts  to  bridge  the  gap  between  the  mathematical  theory  of 
stochastic  processes  and  the  numerical  results  useful  to  the  manage- 
ment science  practitioner. 

Waiting  lines  and  congestion  are  common  problems  encountered  in 
almost  everyone's  daily  life.  A queue  can  develop  in  any  service  system 
whenever  the  immediate  demand  exceeds  the  system's  capacity.  One  of  the 
problems  in  designing  or  controlling  such  systems  is  achieving  the 
proper  balance  between  the  level  of  service  and  the  amount  of  waiting 
which  might  occur.  Formal  analysis  may  be  appropriate  if  the  service 
involves  very  expensive  equipment,  if  the  costs  of  waiting  are  extremely 
high,  or  if  the  decision  involves  many  similar  operating  systems  where 
the  combined  costs  of  service  or  waiting  could  be  excessive. 

The  mathematical  theory  of  queues  can  often  provide  some  insight 
into  the  behavior  of  an  actual  or  proposed  operating  system.  Ideally, 


This  notation  is  the  standard  Kendall  designation  for  queueing  systems: 
Markovian  (or  Poisson)  arrivals/General  service  time  distribution/^ 
(single)  server. 
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the  designer  would  prefer  a mathematical  model  or  algorithm  that  would 
determine  the  optimal  system  design,  given  the  demand  or  arrival  pattern, 
the  costs  of  various  service  levels,  and  the  costs  associated  with 
customer  waiting  or  delay  times.  Unfortunately,  there  is  no  general 
optimization  technique  for  queueing  models.  However,  there  are  numerous 
predictive  models  that  allow  an  analyst  to  determine  some  operating 
characteristics  of  a specific  system.  Prospective  levels  of  service 
can  be  analyzed  one  at  a time  and  an  evaluative  model  can  determine  the 
total  cost  of  service  and  waiting  for  each  alternative. 

The  mathematical  formulas  for  such  operating  characteristics 
that  are  often  available  to  the  practitioner  apply  only  to  queueing 
systems  that  are  in  statistical  equilibrium.  These  steady  state  results 
are  appropriate  for  a system  where  the  server  can,  on  the  average, 
serve  customers  faster  than  they  arrive  and  where  sufficient  time  has 
passed  to  cancel  the  effects  of  the  system's  initial  condition.  If 
customers  arrive  at  the  system  faster  than  the  server  can  handle  them, 
or  if  the  behavior  of  the  system  must  be  analyzed  for  the  start-up 
period,  then  the  time-dependent  (i.e.,  transient)  operating  character- 
istics are  important.  Lee  [1966,  p.  26],  in  an  opinion  probably  shared 
by  many  practitioners,  has  stated  "...  the  first  working  rule  of  queueing 
theory:  time-dependent  solutions  to  queueing-models  are  either  unobtain- 
able or  unmanageable." 

Transient  solutions  that  are  available  in  the  queueing  literature 
are  usually  expressed  in  terms  of  the  Laplace  transform  of  the  operating 
characteristic.  Bhat  [1969,  p.  B284]  reiterates  the  problem  and  mentions 

a possible  solution: 

•"  • ' ■ - 
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"Plainly  speaking,  the  results,  given  in  terms  of  trans- 
forms, very  often  with  more  than  one  argument,  fail  to  make 
sense  to  an  applied  researcher.  Numerical  inversion  of 
transforms  ...  is  an  answer  to  this  problem.  But  at  this 
stage,  the  inversion  methods  are  either  not  sophisticated 
enough  to  handle  the  more  complex  situations  or  do  not  appeal 
to  the  applied  researcher." 

The  present  research  uses  an  accurate  transform  inversion  technique  to 
produce  numerical  results  for  transient  expected  system  load.  A 
practitioner  can  use  these  tables  to  analyze  a wide  range  of  M/G/l 
systems  by  referring  directly  to  Chapter  4.  Chapter  2 discusses  the 
transform  relationships  from  which  we  start,  and  Chapter  3 discusses 
the  inversion  technique  to  be  used. 


1.2.  Methodology 

The  operating  characteristic  studied  in  this  disseratlon  is  time- 
dependent  server  load.  The  server  load  at  epoch  t,  denoted  W(t),  is 
equal  to  the  sum  of  the  service  times  of  customers  waiting  in  the 
queue  plus  any  remaining  service  time  of  a customer  being  served.  The 
quantity  W(t)  is  also  called  virtual  waiting  time,  because  it  is  the 
time  that  a hypothetical  customer  arriving  at  epoch  t would  have  to 
wait  before  beginning  service.  Our  objective  is  to  tabulate  the 
expected  value  of  server  laod,  E[W(t)],  for  various  epochs  t and 
various  system  parameters  (initial  load,  arrival  rate,  and  service  time 
distribution) . 

Our  study  of  M/G/l  queueing  systems  can  be  embedded  in  the  follow- 
ing general  framework.  Let  {X(t),  t > 0}  be  a stochastic  process 


t 
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with  stationary  independent  increments  (a  Levy  process)  whose  sample 
paths  have  no  negative  jumps.  Let  W(t)  be  this  same  process  modified 
by  a reflecting  barrier  at  zero.  If  X(t)  = S(t)  - t,  where 
{S(t),  t > 0}  is  a compound  Poisson  process  with  positive  jumps,  then 
W(t)  is  the  server  load  process  for  an  M/G/l  queue.  (The  jumps  of 
S(t)  occur  at  customer  arrival  epochs  and  the  jump  sizes  are  service 
times.)  We  shall  discuss  two  other  choices  for  the  X process  which 
lead  to  W processes  that  provide  bounds  or  approximations  for  the 
M/G/l  server  load  process.  In  one  instance,  we  take  X(t)  to  be 
Brownian  motion,  and,  in  the  other,  we  take  X(t)  “ G(t)  - t,  where 
G(t)  is  a gamma  process  (a  Levy  process  whose  increments  are  gamma 
distributed) . 

In  Chapter  2 we  present  a Laplace  transform  result  for  E[W(t)] 
that  covers  both  these  two  cases  and  the  queueing  (compound  Poisson) 
case.  In  terms  of  its  application  to  queueing  processes,  this  result 
is  valid  for  a system  that  begins  operation  either  with  or  without  an 
initial  backlog  of  work  and  that  has  an  average  arrival  rate  less  than, 
equal  to,  or  greater  than  its  average  service  rate. 

There  is  no  general  closed-form  expression  for  the  exact  inverse 
of  the  Laplace  transform  of  mean  server  load,  but  numerical  methods  can 
be  used  to  obtain  an  approximation  of  E[W(t}]  at  a specific  epoch  t. 
The  numerical  technique  employed  in  this  research  computes  E[W(t)]  by 
taking  a linear  combination  of  the  Laplace  transform  function  evaluated 
at  appropriate  values  of  its  argument.  Theoretically,  a more  accurate 
approximation  of  E[W(t)]  can  be  obtained  by  using  more  evaluations  of 
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the  Laplace  transform  expression,  but,  when  using  an  electronic  computer 
for  the  computations,  its  finite  word  length  places  a limit  on  the 
accuracy  that  can  be  achieved  in  the  E[W(t)]  approximation.  Investiga- 
tion of  the  technique  using  Laplace  transforms  with  known  inverses 
indicates  that  five  or  six  significant  digits  for  E[W(t)]  can  be 
obtained  efficiently;  this  is  more  than  enough  to  justify  confidence 
in  the  three  or  four  digit  results  shown  in  the  final  E[W(t)]  tables. 

The  computational  procedures  used  in  this  research  could  be  applied 
to  any  M/G/l  queueing  system  whose  service  time  distribution  has  a 
Laplace  transform  that  can  be  evaluated  numerically.  In  this  research 
we  achieve  a balance  between  generality  of  results  and  ease  of  computa- 
tion by  concentrating  on  the  well  known  class  of  M/G/l  queues  whose 
service  times  have  Erlang  (or,  equivalently,  gamma)  distributions.  We 
denote  these  queueing  systems  M/E^/l  and  use  the  Erlang  shape  parameter 
k to  describe  the  service  times.  (Parameter  k is  usually  restricted 
to  positive  integer  values  when  describing  Erlang  distributions,  but 
here  we  allow  k to  assume  any  positive  real  value.)  For  example, 
the  special  case  of  k = 1 corresponds  to  the  exponential  service 
time  distribution  (an  M/M/1  system).  The  present  research  also  examines 
M/Efc/1  queues  with  k less  than  1 and  k greater  than  1,  correspond- 
ing to  service  time  distributions  with  greater  variance  than  the  expo- 
nential and  less  variance,  respectively.  Most  service  time  distributions 
encountered  in  actual  practice  can  be  adequately  approximated  using  the 
very  flexible  Erlang  family. 
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1.3.  Overview  of  the  Subsequent  Chapters 

In  Chapter  2 we  explain  the  sample  path  relationship  between  the 

server  load  process  W(t)  and  the  associated  net  input  process  X(t) 

in  an  M/G/l  queueing  system.  In  this  discussion,  the  net  input  process 

can  actually  be  any  Levy  process  which  has  no  negative  jumps.  Breiman 

[1968]  gives  an  introduction  to  the  theory  of  such  processes,  and 

Blumenthal  and  Getoor  [1968]  provide  a comprehensive  treatment.  We  then 

present  a result  developed  by  Harrison  [1977]  for  the  Laplace  transform 

* 

of  E[W(t)]  when  the  net  input  is  a general  Levy  process.  This 
theoretically  oriented  chapter  finally  examines  a scaling  procedure 
which  facilitates  comparisons  among  the  various  processes. 

Chapter  3 explains  the  Laplace  inversion  technique  developed  by 
Gaver  [1966]  which  gives  approximation  based  on  the  expected  value  of 
an  observational  density  function.  The  accuracy  of  Gaver' s algorithm 
was  improved  by  Stehfest  [1970],  and  we  test  the  technique  using  Laplace 
transform  functions  whose  exact  inverses  are  known.  These  numerical 
results  are  also  compared  with  results  from  Veillon  [1974]  which  were 
obtained  using  other  Laplace  inversion  techniques.  We  then  apply  the 
technique  to  E[W(t)]  in  M/E^/l  queueing  systems,  and  we  compare  our 
results  with  Coleman  [1975]  who  obtained  exact  E[W(t)]  for  the  M/M/1 
case  by  evaluating  sums  of  Bessel  functions.  After  some  additional 
checks  to  ensure  the  accuracy  of  our  approximations,  we  provide  documenta- 
tion of  the  computer  programs  which  generate  the  tables  of  E[W(t)]. 

Gaver  [1966,  1968]  presented  limited  numerical  results  for  transient 
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E [W ( t ) ] in  several  M/G/l  queues,  thereby  demonstrating  the  potential 
usefulness  of  his  technique,  and  Coleman  [1975]  presented  exact  results 
only  for  the  M/M/1  case.  The  major  contribution  of  the  present  research 
is  the  extensive  set  of  tables  in  Appendix  C,  providing  transient  mean 
server  load  in  queues  with  a wide  range  of  traffic  intensities,  initial 
loads,  and  service  time  distributions. 

Chapter  4 explains  how  the  scaled  results  in  the  tables  can  be 
used  by  a practitioner  to  determine  E[W(t)]  in  M/G/l  queues.  Charts 
of  the  Erlang  service  time  distributions  are  presented,  and  simple 
methods  of  interpolation  in  the  tables  are  explained.  Several  sample 
problems  demonstrate  the  entire  procedure.  It  is  intended  that  Chapter  4 
can  be  used  by  a practitioner  without  recourse  to  Chapters  2 or  3. 

The  dissertation  concludes  in  Chapter  5 with  a brief  summary, 
some  qualitative  observations,  and  suggestions  for  further  research. 
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F*(s)  = E(e"sS) 
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e SX  dF(x)  , 


0 < s < 


The  general  approach  used  In  this  research  requires  only  that  this  LST 
can  be  evaluated  numerically.  For  the  specific  cases  to  be  investigated, 
S has  a continuous  distribution  with  density  function  f(*)>  so  the 
LST  reduces  to  the  ordinary  Riemann  integral 


/ e SX  f(x)  dx 


In  particular,  we  examine  M/G/l  systems  with  gamma  distributed 
service  times,  but  we  employ  the  terminology  usually  reserved  for  the 
Erlang  family  of  distributions,  where  the  two  parameters  are  the  mean 
E(S)  and  the  shape  parameter  k.  The  Erlang  density  function  is 


(2.1.1)  f (x)  - x^1  e-kll/E<S) 


x > 0 , 


with  k > 1 restricted  to  integer  values.  However,  here  we  allow  k 
to  assume  any  non-negative  real  value,  requiring  that  the  factorial 
(k— 1) ! in  the  density  function  be  replaced  by  the  gamma  function, 

oo 

T(k)  = / xk  1 e X dx 
0 

Using  the  Erlang  terminology,  the  LST  of  our  service  time  distribution 
is  (see  Drake  [1967,  p.  138]  for  a derivation) 
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F*(s)  - [ 


k + sE(S) 
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k > 0 , 


The  variance  of  an  Erlang  random  variable  is  [ E (S) ] /k.  Thus, 
for  various  Erlang  service  time  distributions  with  the  same  mean,  the 
shape  parameter  k is  inversely  proportional  to  the  variability  of 
those  service  times.  For  example.  Erlang  service  times  with  0 < k < 1 
have  even  more  variability  than  an  exponential  distribution  (k  = 1). 

On  the  other  hand,  for  very  large  values  of  k the  variance  is  very 
small,  and  we  approach  the  case  of  a constant  or  deterministic  service 
time  (an  M/D/1  system)  as  k tends  to  infinity.  Probability  density 
functions  for  0 < k < 1 and  k > 1 are  shown  graphically  in  Figure  4.1. 
Hillier  and  Lieberman  [1974,  p.  417]  state  that  "empirical  service-time 
distributions  can  usually  be  reasonably  approximated  by  an  Erlang 
distribution. " 

The  most  important  descriptive  parameter  for  a queueing  system 
is  its  traffic  intensity  p,  defined  as  the  mean  service  time  divided 
by  the  mean  interarrival  time.  Thus,  for  M/G/l  systems,  we  have 


p = \E(S) 


When  p is  less  than  one,  i.e.,  when  the  average  time  between  arrivals 
is  greater  than  the  average  service  time,  we  say  that  the  system  is 
"stable".  In  this  case,  p is  the  fraction  of  time  that  the  server  is 
busy,  and  it  can  be  interpreted  as  the  system's  "utilization  factor". 
Furthermore,  for  p < 1 the  distribution  of  virtual  waiting  time 


10 


1 


approaches  an  equilibrium  distribution  as  t increases.  Let  W denote 
this  steady-state  waiting  time.  Its  expected  value  is  given  by  the 
Pollaczek-Khintchine  formula, 


I 

[•:  (2.1.3) 


E(W) 


\E(S2) 

" 2(l-p)  ’ 


pel  . 


Our  numerical  results  for  time-dependent  virtual  waiting  time  will 
allow  us  to  observe  how  quickly  this  steady-state  condition  is  approached. 
However,  we  will  also  examine  "unstable"  systems  (p  > 1),  where  the 
queue  length  and  waiting  time  tend  to  increase  without  bound,  so  that 
no  steady-state  conditions  exist. 

As  mentioned  in  the  introductory  chapter,  the  virtual  waiting  time 
or  server  load,  W(t),  represents  the  work  backlog  at  epoch  t,  i.e., 
the  accumulated,  unserviced  demand.  We  define  W(t)  in  terms  of  the 
underlying  work  input  processes,  which  allows  us  to  use  existing  Laplace 
transform  results  for  reflected  Levy  processes.  Sample  path  relationships 
are  shown  graphically  in  Figures  2.1  and  2.2,  illustrating  one  possible 
realization  for  these  stochastic  processes. 

The  input  process  S(t)  represents  the  amount  of  customer  work 
that  arrives  during  the  interval  [0,t]  and  is  defined  by 


S(t)  - Sj  + •••  + sA(t)  , t » 0 . 


This  compound  Poisson  process  has  mean  E[S(t)]  = pt  and  variance 
Var [S(t) ] = \E(S2)t. 
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The  net  input  process,  X(t),  is  defined  by 


X(t)  * S(t)  - t , t > 0 . 

If  the  first  customer  arrives  at  t * 0,  then  X(t),  sometimes  called 
pseudo-server  load,  is  identical  to  the  server  load  process  as  long  as 
the  server  remains  busy.  The  similarity  ends  when  the  server  first 
becomes  idle.  A general  expression  for  W(t)  in  terms  of  X(t) 
requires  that  we  account  for  server  idleness. 

Let  I(t)  represent  the  cumulative  idleness  of  the  server  during 
the  interval  [0,t],  Since  the  server  works  at  a unit  rate,  the 
potential  work  output  is  t units  during  the  same  interval,  and  the 
actual  work  output  is  t - I(t).  Assuming  no  initial  work  load,  i.e., 

W(0)  = 0,  the  available  work  during  the  interval  [0,t]  is  S(t), 
so  the  remaining  work  load  W(t)  can  be  expressed  as  the  difference 
between  the  work  input  and  the  actual  work  output,  i.e.,  W(t)  * S (t)-[ t-I (t) ] , 
or  W(t)  = X(t)  + I(t).  The  sample  path  relationships  can  be  seen  in 
Figure  2.1.  When  X(t)  is  equal  to  its  infimum,  the  server  is  idle. 

When  X(t)  is  greater  than  its  infimum  (because  some  work  has  arrived), 
the  server  is  busy  serving  customers.  The  infimum  becomes  more  negative 
only  when  the  server  is  idle.  Therefore,  the  cumulative  idleness  can 
be  expressed  as 

I(t)  - -inf[X(T)  : 0 < t < t]  , t > 0,  if  W(0)  - 0 . 
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Combining  this  with  the  net  input,  we  have  the  server  load  representa- 
tion (no  initial  workload). 


W(t)  - X(t)  - inf [X(t)  : 0 < t < t]  , 


t > 0,  if  W(0)  - 0 


The  same  reasoning  applies  to  the  more  general  case  of  initial 
server  load,  illustrated  in  Figure  2.2,  where  the  underlying  S(t) 
and  X(t)  sample  paths  would  be  identical  to  those  shown  in  Figure  2.1. 
The  initial  server  load  z represents  a specific  amount  of  work  avail- 
able for  service  at  epoch  t ■ 0.  Then  the  available  work  during  an 
interval  [0,t]  is  z + S(t),  and  the  remaining  work  load  W(t)  at 
any  epoch  t is  z + S(t)  - [t  - I (t) ] , or  W(t)  ■ z + X(t)  + I(t). 

The  cumulative  idleness  function  I(t)  is  slightly  more  complicated 
when  W( 0)  > 0.  The  server  is  initially  busy  in  this  case,  so  when 
z + X(t)  is  equal  to  its  infimum,  the  server  is  idle  only  if  that 
inflmum  is  negative.  After  the  initial  busy  period,  the  infimum  of 
z + X(t)  becomes  more  negative  only  when  the  server  is  idle,  and  the 
representation  of  the  cumulative  idleness  function  is  similar  to  the 
case  with  no  initial  load.  Thus,  for  this  general  case  the  server  load 
representation  is 


W(t)  • z + X(t)  - {inffz  + X(t)  : 0 < t < t]}  , 

t > 0,  if  W(0)  ■ z > 0 . 
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I 


i 


* 

T j 

Our  main  objective  is  to  compute  the  expected  value  of  W(t), 
where  the  expectation  is  taken  with  respect  to  a probability  distribu- 
tion over  all  possible  sample  paths.  We  adopt  the  notation  Ez [W ( t ) ] 

and  E [I (t)  ] for  expected  server  load  and  expected  cumulative  idleness, 
z 

respectively,  conditional  on  initial  work  load  W(0)  = z.  From  our 
discussion  of  the  sample  path  relationships  it  follows  that 

E [W(t)]  = E [z  + S(t)  - t + I(t)]  , 

z z 

(2.1.3) 

E [W(t) ] “z+pt-t+E  [I(t) ) . 
z z 

Thus,  in  M/G/l  queueing  systems  the  mean  server  load  can  be  separated 
into  four  additive  terms:  initial  work  load,  new  work  input,  potential 
work  output,  and  cumulative  idleness.  Another  useful  observation  is 
that  E^ [ I ( t ) ] must  be  zero  in  the  interval  from  t * 0 to  t * z, 
i.e.,  it  is  impossible  for  the  server  to  become  idle  before  the  initial 
work  load  has  been  serviced.  In  Chapter  3 this  property  is  used  to 
provide  a check  on  the  accuracy  of  our  numerical  results. 

Thus  far  the  sample  path  relationships  and  expectations  have 
been  discussed  in  the  context  of  M/G/l  queueing  systems.  However,  we 
also  calculate  Ez[W(t)]  for  processes  where  the  same  relationships 
apply,  but  where  S(t)  is  not  compound  Poisson.  These  related  processes 
provide  bounds  or  approximations  for  M/G/l  mean  server  load. 
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2.2.  The  Laplace  Transform  Result  for  Expected  Server  Load 

All  of  the  stochastic  processes  for  which  numerical  values  are 
calculated  here  can  be  discussed  simultaneously  in  the  following 
unified  framework.  Let  X ■ {X(t),  t > 0}  be  a process  with  stationary. 
Independent  increments  (an  infinitely  divisible  or  Levy  process)  and  no 
negative  jumps,  and  define 


-E[X(t)]  - nt  , 


where  -»  < p.  < - , 


(2.2.1) 


Var[X(t) ] -at. 


where  0<  a < - , 


for  t > 0.  According  to  the  standard  results  for  Levy  processes,  the 
Laplace  transform  of  X(t)  has  the  form 


E[«-"(t)]  - e-*<8)t 


for  s > 0 and  t > 0 , 


and  the  exponent  function  $(•)  is  convex  with 


(2.2.2) 


$(0)  - 0,  $'(0)  » p.  , and  $"(0)  - o , 


cf.  Harrison  [1977]  and  Takacs  [1967].  It  can  be  shown  that  for  each 
s > 0 there  exists  a unique  co(s)  > 0 such  that 


(2.2.3) 


4>[co(s>  ] - s 


J 
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In  Che  previous  section  we  discussed  Che  sample  path  relationships 
where  process  W is  obtained  from  X by  imposing  a reflecting  barrier 


■I 


at  zero.  Using  the  notation  Ez[W(t)]  = E[W(t) |W(0)  = z],  let  Pw(z,s) 

denote  the  Laplace  transform  of  E [W(t>],  that  is, 

z 

oo 

Pu(z,s)  - / e'St  E [W(t) ] dt  , s > 0 . 

w o z 

Harrison  [1977]  developed  a simple  formula  for  the  Laplace  transform  of 

E [W(t)],  where  p.  is  unrestricted  in  sign: 
z 


(2.2.4) 


Pw(z,s) 


-w(s)z 

e 

sw(s) 


In  Chapter  3 we  use  this  formula  to  obtain  accurate  approximations 

of  E [W(t) ] for  specific  epochs  t,  initial  loads  z,  and  various 
z 

net  input  processes  X.  To  achieve  the  desired  accuracy  the  numerical 
Inversion  technique  requires  that  Pw(z,s)  be  computed  for  34  values 
of  s in  order  to  obtain  Ez(W(t)]  at  a single  epoch,  and  each  evalua- 
tion of  P^(z,s)  requires  that  the  functional  equation  (2.2.3)  be 
solved  to  obtain  co(s). 

The  ease  with  which  co(s)  can  be  calculated  depends  upon  the  form 
of  the  exponent  function,  and  the  exact  form  of  4>(*)  depends  upon 
the  process  X.  If  we  specify  X(t)  ■ S(t)  - t,  where  S(t)  is  a 
compound  Poisson  process,  then  it  can  be  shown  that 


$(s)  - s - X[1  - F*(s) ] , 
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In  the  previous  section  we  discussed  the  sample  path  relationships 
where  process  W is  obtained  from  X by  imposing  a reflecting  barrier 
at  zero.  Using  the  notation  E^tVKt)]  = E[W(t)|w(0)  = z],  let  Pw(z,s) 
denote  the  Laplace  transform  of  E [W(t)],  that  is. 


Pu(z,s)  = / e"&t  E [W(t) ] dt  , 
W 0 z 


s > 0 


Harrison  [1977]  developed  a simple  formula  for  the  Laplace  transform  of 

E [W(t)],  where  p is  unrestricted  in  sign: 
z 


(2.2.4) 


Pw(z,s) 


L _ J±_  + 
s 2 


a-u(s)z 

sw(s) 


In  Chapter  3 we  use  this  formula  to  obtain  accurate  approximations 

of  Ez[W(t)]  for  specific  epochs  t,  initial  loads  z,  and  various 

net  input  processes  X.  To  achieve  the  desired  accuracy  the  numerical 

inversion  technique  requires  that  P (z,s)  be  computed  for  34  values 

w 

of  s in  order  to  obtain  Ez[W(t)]  at  a single  epoch,  and  each  evalua- 
tion of  Pw(z,s)  requires  that  the  functional  equation  (2.2.3)  be 
solved  to  obtain  oo(s). 

The  ease  with  which  w(s)  can  be  calculated  depends  upon  the  form 
of  the  exponent  function,  and  the  exact  form  of  <*>(•)  depends  upon 
the  process  X.  If  we  specify  X ( t ) * S(t)  - t,  where  S(t)  is  a 
compound  Poisson  process,  then  it  can  be  shown  that 


4>(s)  - s - \[1  - F*(s)  ] , 


cf.  Prabhu  [1965,  p.  70]  and  Takacs  [1967,  p.  59].  Using  the  terminology 

of  M/G/l  queueing  theory  X is  the  average  arrival  rate  and  F*  C * ) is 

the  Laplace  transform  of  the  service  time  distribution.  Since  S(t) 

2 

has  mean  XE(S)t  and  variance  XE(S  )t,  it  follows  that  the  parameters 

2 2 

of  X defined  by  equations  (2.2.1)  are  = 1-p  and  o * \E(S  ).  If 
F*(#)  can  be  evaluated  numerically,  then  the  properties  (2.2.2)  of 
#(•)  indicate  that  the  functional  equation  $[oo(s)]  = s can  be  solved 
efficiently  using  an  elementary  one-dimensional  search  technique.  In 
Chapter  3 we  discuss  our  use  of  the  Newton-Raphson  method  to  obtain 
<o( s)  for  M/E^/l  queueing  systems. 

For  the  special  case  of  an  M/M/1  system  the  LST  of  the  service 
time  distribution  is  obtained  by  setting  k **  1 in  equation  (2.1.2), 
i.e.,  F*(s)  * 1/ [1  + sE(S) ] . For  a specified  value  of  s the  functional 
equation  <i>[a>(s)]  * s is  a quadratic  function  of  a>(s) , and  the  positive 
solution  is 


(2.2.5) 


-U  - x[s+e(s)  i ) + m - x[s+E(s; 

2E(S) 


+ 4sE(S) 


The  M/M/1  system  is  the  only  queue  studied  here  that  has  an  analytic 
solution  to  (2.2.3).  In  general  the  M/D/1  and  M/E^/l  cases  will 
require  a search  to  determine  co(s) . The  formulas  for  these  cases 
are  summarized  in  Table  2.3  below. 

Two  other  choices  for  X(t)  yield  reflected  processes  W(t) 
which  provide  bounds  or  approximations  for  M/G/l  server  load.  The 
first  case  is  the  Wiener  process,  which  has  been  used  as  a model  for 


a variety  of  physical  processes  since  its  original  development  as  a 


Gamma  M/G/l  Queueing  Systems  Brownian 
Input  Motion 
Process  Process 


TABLE  2.3.  Summary  of  Characteristics  for  the  Processes  under  Study. 


model  for  Brownian  motion.  Gaver  [1968]  proposed  that  this  diffusion 
process  could  be  used  as  an  approximation  for  the  net  Input  process 
of  an  M/G/l  queue  by  equating  the  first  and  second  moments  of  the  two 


processes. 

The  Brownian  motion  process  is  denoted  by 
X^t)  = o£(t)  - nt  , 

where  £(t)  is  a Wiener  rocess  whose  stationary,  independent  increments 

have  a normal  distribution  with  mean  zero  and  unit  variance.  It  follows 

2 

that  X_ (t)  has  mean  -|it  and  variance  at.  As  previously  noted, 

the  net  input  process  of  an  M/G/l  queue  has  mean  E[X(t)]  = (p-l)t  and 

2 

variance  Var[X(t)]  = \E(S  )t.  Thus,  using  Xg(t)  to  approximate  X(t) 

2 2 

we  would  take  -p.  = p-1  and  a = \E(S  ).  The  exponent  function  for 

1 2 2 

Brownian  motion  is  $(s)  = (is  + y a s , cf.  Takacs  [1967,  p.  81];  thus 
the  solution  of  $[<a(s)]  * s can  be  calculated  using  the  quadratic 
formula. 

The  second  non-queueing  input  process  that  we  consider  as  an 
approximation  is  a gamma  input  process,  denoted 

XG(t)  - G(t)  - t , 

where  G(t)  is  a process  whose  stationary,  independent  increments  have 
a gamma  distribution.  The  gamma  density  function  is 


m 


2 

so  G(t)  has  mean  (a/p)t  and  variance  (a/p  )t.  It  follows  directly 

2 

that  X (t)  has  mean  (a/p  - l)t  and  variance  (a/p  )t.  Thus,  using 

Vj 

Xg(t)  to  approximate  the  net  input  process  of  an  M/G/l  queue  we  would 

2 2 

choose  a and  (3  such  that  a/p  - 1 = p-1  and  a/p  = \E(S  ).  The 
exponent  function  for  the  gamma  input  process  is  <f>(s)  = s - a log(l  + s/p), 
cf.  Takacs  [1967,  p.  66];  the  solution  of  the  functional  equation  (2.2.3) 
requires  a search  technique. 


2.3.  Scaling  the  Reflected  Levy  Processes 

In  this  section  we  present  a method  for  normalizing  the  processes 
of  interest  by  using  a simple  linear  tranformation  for  both  the  epochs 
and  the  server  load.  This  technique  allows  us  to  reduce  the  tabulations 
required  to  describe  actual  operating  systems  and  also  facilitates 
comparisons  among  different  processes  by  adopting  a common,  normalized 
time  scale. 

Consider  an  M/G/l  queueing  system  (as  originally  described  in 

Section  2.1)  with  service  times  S^,  S2,  ...  having  distribution 

2 

function  F(-)  with  mean  E(S)  and  second  moment  E(S  ).  Let 
{A(t);  t > 0}  be  the  Poisson  arrival  process  with  mean  arrival  rate  X. 
The  net  input  process  is 

X(t)-s1  + ... +sA(t)-«  , 
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and  the  server  load  process  W(t)  is  obtained  from  [z  + X(t)]  by 
the  reflection  mapping  described  in  Section  2.1,  where  z is  the 
initial  server  load  W(0). 


We  define  a new  scaled  process  in  terms  of  the  original  queueing 


process 


X (t)  = aX(bt) 


This  scaled  net  input  process  has  the  form 


(2.3.1)  X*(t)  = S*  +•••  + S*  - ct  , 


t > 0 , 


* * 

where  c ■ ab,  A (t)  = A(bt),  and  * aS^  Note  that  the  random 
* * 

variables  have  distribution  function  G (x)  ■ F(x/a),  and  that 

* * * 
A (t)  is  a Poisson  process  with  mean  arrival  rate  X * bX.  Let  p. 

2 

and  a * be  defined  by 


E[X*(t)]  « -A  , 


(2.3.2) 


Var  [X*(t)  ] = a2,  t 


We  accomplish  our  normalization  by  choosing  a and  b so  that 

i*i  2 

| M-  I * 1 and  * 1. 

By  substituting  aX(bt)  for  X*(t)  on  the  left  hand  side  of 
equations  (2.3.2)  and  then  using  the  parameters  of  process  X as 
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1 


m 

f 

£ 


scallne  «*«  ’ u u e*’y  * shov  th 


and 


Ifwe  deflne  w*  ^ 

SCal6d  *««al  se  38  the  reflectlon  ^ 

that  **,..*  Ioad  is  z*  , +X*(t)] 

f ^ * a»0»t;.  az»  then  lt  , ’ here  the 

Process  < That  is,  th  ±s  easy  to  , 

8 ld6n«c ai  ^ th  refleCtl°n  of  th 

Process.  If  the  scaled  Ve  he  soaled  net 

Process,  th  6 l6t  ^ represent  a ^ ^ **  °ri8inal  Serv 

the  scaled 


£z*fv*(t 2 

o * t*^ 

4^  J * 


*"■  “e  - scaled  “ ' 2 ' ^ ■*  . 

,Ueu«»«  syste„  “ «rver  „ 

^ d<«u  j„7rs  the  *n  °ri8i"*‘ 

rtas  — -s  r scaiM  ~ 
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that  ,a  00<i  hueue, „„  Tke  craff, 
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X&(S) 


’ " 'fi 


litfjft'-flfudl  Mill 
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■ ,w , 


but  Che  variance  parameter  Is  different: 


a'2  - X*  E(S ' 2)  = X'  (pr)2  E(S2)  * YT  ^(S2)  = TT  o2 


This  second  system  and  the  original  are  related  by 


x’(t)  -^(yt)  , 


or  equivalently  by 


X'(t) 


E(S')  y,E(S) 
E(S)  ME(S’) 


t) 


In  the  context  of  M/E^/l  queueing  systems,  the  second  queueing  system 

has  the  same  p and  same  Erlang  shape  parameter  k as  the  original, 

2 2 

but  the  different  values  for  a and  ct'  require  simple  transforma- 
tion of  the  epoch  scale  and  waiting  time  scale. 

Returning  to  the  scaled  process  which  was  defined  in  terms  of 
the  original  queueing  process,  we  now  express  the  original  process  in 
terms  of  the  second  process: 


X*(t) 


- M */2 i 


X(^  t) 


Mi 

2 X' 
a 


X'  (T  *2  t} 


CT*  |JL 


Once  p (or  |J.)  has  been  specified  and  a specific  form  of  the  service 
time  distribution  (i.e.,  a specific  Erlang  shape  parameter  k)  has  been 
determined,  then  we  are  free  to  choose  the  variance  parameter  (i.e., 
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the  time  scales)  for  the  queueing  process  from  which  X*  will  be 
obtained. 

This  research  tabulates  EzA[W*(t*)]  for  various  values  of  p, 
k,  and  z*.  and  one  can  obtain  mean  server  load  for  an  actual  queueing 
system  by  selecting  the  table  with  the  closest  (p,  k,  z*)  combination 
or  by  interpolating  between  two  tabulated  scaled  processes,  and  then 
applying  the  transformation: 


Ez[W(t)] 


Vf 


a 


t)] 


Since  p.  - 1-P,  this  particular  normalization  cannot  be  used  for  systems 

with  p * 1,  i.e.,  p.  = 0.  In  this  case  we  tabulate  mean  server  load 

for  Erlang  queueing  systems  with  X and  E(S)  chosen  so  that  the 

2 2 

variance  parameter  a m XE(S  ) * 1.  Let  superscript  zero  indicate 

2 

parameters  for  the  tabulated  system.  Then  a = 1 requires  that 


.0  k+1 

K “ T 


and 


E(S°) 


k 

k+1  • 


Mean  server  load  for  an  actual  operating  system  with  shape  parameter  k 
and  mean  arrival  rate  X can  be  obtained  from  the  tabulated  values  as 
follows : 


EztW(t)] 


X 


! 0[“0<75  ‘>1 

z X 


k+1 

kX 


E 0tW  ^k+i  ’ 
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Ijl 
. 


- I 


!!  I 


I 


The  tabulated  values  for  the  p * 1 case  correspond  to  an  actual 
queueing  system,  whereas  the  tabulated  values  for  the  p f1  1 cases 
represent  a stochastic  process  with  "potential  work  output  rate"  of 
c * ab  ■ 1/  | M- 1 » Interpreted  from  equation  (2.3.1).  For  the  p<  1 
cases,  mean  server  load  In  queueing  systems  approaches  the  Pollaczek- 
Kh in t chine  values,  equation  (2.1.3).  Thus,  the  tabulated  values  for 
all  scaled  processes  approach 


M . 


2TmT 


thereby  allowing  comparisons  concerning  approach  to  equilibrium. 
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CHAPTER  3 


r 


NUMERICAL  INVERSION  OF  THE  LAPLACE  TRANSFORM 

, 

3.1.  Approximation  by  an  Expected  Value 

In  the  previous  chapter  we  cited  an  expression  for  the  Laplace 
transform  for  the  three  processes  of  interest.  For  the  two  special 

I cases  of  Brownian  motion  and  the  M/M/1  queue,  analytic  methods  can 

be  applied  and  exact  numerical  solutions  can  be  obtained.  However, 
in  general  it  is  necessary  to  use  a method  for  numerical  inversion  of 
the  Laplace  transform  in  order  to  evaluate  the  process  of  interest. 

In  this  chapter  we  describe  such  a method  and  its  implementation. 

The  first  two  sections  of  the  chapter  describe  the  method  lor 
approximate  transform  inversion  using  an  expected  value.  In  the  third 
section  we  present  numerical  results  for  test  cases  where  this  method 
is  applied  to  Laplace  transforms  whose  exact  inverses  are  known,  and 
we  describe  the  computer  routines  that  apply  the  inversJDn  method  to 

I prepare  tables  of  E[W(t)]. 

The  method  used  in  this  research  for  numerical  inversion  of  the 
Laplace  transform  was  originally  developed  by  Gaver  [1966].  Although 
there  are  other  techniques  which  could  have  been  implemented,  this 
particular  method  was  chosen  because  it  had  been  applied  successfully 
to  the  Laplace  transform  expression  for  E[W(t)]  in  the  M/M/1  and 
M/G/l  queues  by  Gaver  [1966,  1968]. 

The  problem  of  inverting  the  Laplace  transform  can  be  summarized 
as  follows.  We  have  a function  of  Interest,  P(t),  for  which  no  closed- 
form  analytic  expression  is  known,  so  that  it  cannot  be  evaluated 
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numerically.  However,  we  do  know  the  expression  for  its  Laplace 
transform, 


p (s)  = / e st  P(t)  dt 


Specifically,  in  this  research  the  function  or  process  of  interest  is 
E[W(t)]  and  the  Laplace  transform  expression  was  cited  in  Chapter  2. 

In  general,  we  wish  to  tabulate  P(t)  for  various  values  of  t. 

The  method  presented  here  computes  an  approximate  value  of  the 
function  at  a specified  point  t'.  Theoretically  this  approximation 
can  be  computed  as  precisely  as  desired.  Consider  observing  the  function 
of  interest  at  a random  time  T that  has  density  function  f(t) 
such  that  the  values  of  T are  concentrated  near  t',  as  shown  in 
Figure  3.1.  Then  P(t'),  an  approximation  of  P(t'),  can  be  expressed 


(3.1.1) 


P(t’)  - / P(t)  f(t)  dt 

0 


Now  the  problem  is  one  of  finding  an  observational  density  function 
f(t)  so  that  the  right  hand  side  of  (3.1.1)  can  be  expressed  in  terms 
of  p(s). 

Gaver  [1966]  examined  such  a family  of  density  functions  that 
are  well-suited  to  numerical  computation.  Let  fn(t;  a)  be  t*ie  density 
function  for  random  variable  T,  with  parameters  n and  a,  as  follows: 


0.1.2)  f„(t,  .)  - « ^T)T  » - e~‘t>n  e'"at  • *»0.  .-1,2... 

1 I 

Gaver  [1966]  showed  that  f (t;  a)  has  the  following  properties: 


modal  value  ~ — In  2 , 

~ a 

„ 1 2n¥2 

Var(T)  ~ - (2n-l)  (4n+l)  * 

O 

For  Increasing  values  of  n,  f (t;  a)  becomes  more  sharply  peaked  at 

n 

t - 1/a  In  2.  Using  f (t;  a)  as  the  observational  density  function 


in  (3.1.1),  we  let  the  approximation  of  P(t')  be 

oo 

(3.1.3)  P - / P(t)  f (t;  a)  dt  , 

n 0 n 

where  parameter  a is  chosen  such  that  a * 1/t'  In  2,  thereby 
concentrating  the  values  of  random  variable  T near  t'  as  desired. 

By  substituting  the  observational  density  (3.1.2)  in  the  approximating 

£ fj 

expression  (3.1.3)  and  then  expanding  (1  - e ) by  the  binomial  theorem, 
we  obtain 

(3.1.4)  P - a 2 (")  (-1)1  p[(n+i)a]  . 

n n!  (n-1) ! issQ  i , 

Thus,  P^,  an  approximation  of  P(t'),  is  based  on  a linear  combination 
of  the  Laplace  transform  p(s)  evaluated  at  n+1  different  values 
of  s.  Theoretically,  the  sequence  {P^;  n = 1,  2,  3,  ...}  converges 
to  P(l/a  In  2),  i.e.,  to  P(t').  For  any  finite  n,  we  can  calculate 
P using  (3.1.4),  and  because  Var(T)  -►  0 as  n -*■  00  we  can  obtain 

as  good  an  approximation  as  desired  by  using  a large  enough  value  of  n. 
Practically,  p(s),  the  Laplace  transform  of  P(t),  must  be  evaluated 
at  n+1  values  of  s in  order  to  evaluate  P^;  in  some  applications 
much  time  may  be  required  to  make  each  evaluation  of  the  Laplace  transform 
expression  and  so  the  cost  of  obtaining  the  desired  precision  in  Pr  may 
be  prohibitive.  Also,  for  large  n the  factorial  terms  in  (3.1.4)  are 
integers  with  too  many  digits  to  be  precisely  evaluated  on  computers 
with  finite  word  length,  and  operations  involving  these  large  numbers 


3 


can  result  in  considerable  rounding  errors.  From  the  standpoint  of 

computational  efficiency  it  is  important  to  improve  the  accuracy  of 

the  approximate  inverse  without  requiring  evaluation  of  P for 

n 

extremely  large  values  of  n. 


3.2.  Improving  the  Accuracy  of  the  Approximation 

Gaver  [1966]  showed  that  each  can  be  represented  as  an 

asymptotic  expansion,  i.e., 

“l  a2  a3 

(3.2.1)  P « P(t')  + — + + -£  + ••• 

n n 4 J 

n n 

where  the  error  components  depend  only  on  t'  and  not  on  n.  An 

improved  approximation  of  P(t')  can  be  obtained  by  taking  a linear 
combination  of  a set  of  the  P approximations,  where  the  values  of 
n are  integer  powers  of  2.  This  method,  extrapolation  to  the  limit, 
ensures  that  some  of  the  error  terms  in  (3.2.1)  are  cancelled  out, 
cf.  Gaver  [1966]  and  Henrici  [1964]. 

Stehfest  [1970]  demonstrated  an  improved  calculation  method 
that  uses  a linear  combination  of  an  even  number  of  the  P^  approxima- 
tions; this  method  cancels  even  more  of  the  error  terms  than  Gaver’ s 
method  (extrapolation  to  the  limit).  Stehfest  combined  the  coefficients 
for  the  linear  combination  of  P^  with  the  coefficients  of  p(s)  in 
(3.1.4)  so  that  Fa,  the  approximation  of  P(t'),  could  be  expressed 
directly  as  a linear  combination  of  p(s): 
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(3.2.2) 


In  2 „ „ ,ln  2 

—rr  2 V p(— r-  i) 

c i=l 


Here  N must  be  even  and  the  combined  coefficients  are 


(3.2.3) 


(N/2)+i 


Min(i,-|) 


kN/2  (2k)! 


k[l+l]  (f-k)!  k!  (k-1) ! (i-k)!  (2k-l) ! 


For  a theoretical  comparison  of  the  original  Pr,  extrapolation 
to  the  limit,  and  Fa,  consider  the  three  methods  where  the  number  of 
values  of  p(s)  is  32  in  each  case.  Thirty-two  evaluations  of  the 
p(s)  could  be  used  to  compute  P^,  and  from  (3.2.1)  the  first  error 
term  would  be  a /31.  On  the  other  hand  the  thirty-two  values  of 
p(s)  could  be  used  to  compute  P , P2>  P^,  Pg,  and  P^^,  and  it 
could  be  shown  that  extrapolation  to  the  limit  would  yield  an 
approximation  of  P(t')  where  the  first  four  error  terms  of  equa- 
tion (3.2.1)  would  cancel  completely.  The  resulting  approximation 


would  be  P(t')  + a^/2  + 


Finally  using  Stehfest's  method 


with  N * 32,  the  first  fifteenerror  terms  in  equation  (3.2.1) 
would  cancel,  and  it  could  be  shown  that  the  approximation  would 
be  Fa  = P(t')  - a^/16!  + •••  . This  theoretical  superiority  of 
Stehfest's  method  was  verified  numerically  in  preliminary  investi- 
gations using  functions  with  known  inverses. 
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3.3.  Computer  Implementation  of  the  Technique 


In  this  research,  the  algorithm  LINV  developed  by  Stehfest 
[1970]  was  translated  from  ALGOL  to  FORTRAN  IV.  Some  changes  were 
made  in  evaluating  the  V of  equation  (3.2.3)  to  take  advantage  of 
cancelling  the  factorial  terms  and  thereby  avoid  rounding  errors.  Some 
preliminary  numerical  investigations  used  BASIC  and  single  and  double 
precision  FORTRAN  IV.  The  computer  programs  listed  in  Appendix  B 
specify  double  precision  for  all  non- integer  variables,  and  the  "AUTODBL" 
option  available  with  the  IBM  FORTRAN  H-Extended  Compiler  was  used  in 
all  of  the  research  reported  here.  With  this  extended  precision  there 
are  approximately  36  significant  decimal  digits  for  the  non-integer 
variables. 

Stehfest  [1970]  published  a table  of  results  applying  LINV  to 
six  transforms  whose  inverses  are  known,  using  8 digit  arithmetic  and 
n = 10.  Figure  3.2  shows  che  exact  values  of  the  six  functions,  the 
approximations  from  this  research  using  N = 34  with  36  digit  arithmetic, 
and  the  original  Stehfest  approximations  using  N = 10.  We  observe  that 
there  are  at  least  six  correct  significant  digits  in  our  approximations 
(using  N = 34)  for  each  of  the  test  functions. 

Since  the  LINV  approximation  (Fa)  is  based  upon  the  original  P , 
it  should  also  be  more  accurate  when  more  values  of  the  Laplace  transform 
p(s)  are  used,  i.e.,  if  N is  very  large.  However,  for  large  values 
of  N the  rounding  errors  in  evaluating  V^  can  affect  the  accuracy 
of  the  results.  This  problem  was  investigated  by  applying  LINV  to  the 
six  test  functions  of  Figure  3.2;  approximations  were  obtained  using 
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FIGURE  3.2  INVERSION  OF  TEST  FUNCTIONS 
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N equal  to  30,  32,  34,  and  36.  In  each  case  the  correct  number  of 
decimal  places  was  recorded,  rounding  off  both  the  exact  value  and  the 
approximation  when  making  the  comparison.  Figure  3.3  shows  the 
average  number  of  correct  decimal  places  over  the  ten  values  of  t 
evaluated  in  each  case.  We  observe  that  the  approximations  using 
N = 34  are  at  least  as  good  as  those  using  N = 32  or  N = 36  in 
all  cases.  Thus,  in  the  remainder  of  this  research  we  will  use  N = 34 
in  the  LINV  algorithm.1 

Figure  3.4  compares  the  LINV  algorithm  with  other  inversion 

techniques.  The  top  section  of  this  table  reproduces  some  results  by 

2 

Dubner  and  by  Veillon  (1974)  ; Dubner's  approximations  are  based  on 
500  values  of  p(s),  while  Veillon's  method  uses  64  values.  We  observe 
that  in  all  cases  the  LINV  algorithm  using  34  values  of  p(s)  is  at 
least  as  accurate  as  their  techniques.  The  bottom  section  of  Figure  3.4 
compares  approximations  of  a second  function;  we  observe  that  once  again 
the  LINV  results  are  as  accurate  as  the  results  obtained  by  other 
techniques. 

Figures  3.2,  3.3  and  3.4  indicate  that  LINV  can  provide  very 
accurate  approximations  for  the  test  functions,  but  it  is  also  important 
to  investigate  LINV's  accuracy  using  the  function  of  interest  in  this 

1In  preliminary  investigations  using  FORTRAN  IV  double  precision 
(approximately  17  digit  arithmetic)  without  the  AUTODBL  option,  the 
most  accurate  approximations  were  obtained  using  N * 18. 

2 

Reprinting  privileges  for  the  Stehfest  results  (Comm,  of  the  ACM,  Vol.  13, 
No.  1,  Copyright  1970)  and  the  results  published  by  Veillon  (Comm,  of 
the  ACM,  Vol.  17,  No.  10,  Copyright  1974)  were  granted  by  permission 
of  the  Association  for  Computing  Machinery. 


FIGURE  3.3.  Average  Number  of  Correct  Rounded-off  Decimal 
Places  in  LINV  Approximations  of  Six  Functions 
for  Determination  of  Optimal  Value  of  N 


N,  the  number  of  values  of  the 
Laplace  transform  p(s)  used 
by  LINV  to  estimate  the  function 
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research,  E[W(t)].  Our  earlier  investigations  using  LINV  (not  included 
here)  verified  the  results  for  M/M/1  queues  reported  by  Gaver  [1968], 
but  it  is  more  appropriate  to  compare  our  LINV  results  with  mean 
server  load  obtained  by  some  method  other  than  Laplace  transform  inver- 
sion. Fortunately,  Coleman  [1975]  evaluated  E[W(t)]  in  the  M/M/1 
queue  using  a sum  of  Bessel  functions.  In  Figure  3.5  the  first  three 
columns  show  his  results  at  five  epochs  and  the  LINV  results  using 
equations  (2.2.4)  and  (2.2.5)  for  the  Laplace  transform  expression. 

The  five  significiant  digits  available  in  the  Coleman  results  are 
matched  exactly  when  the  LINV  results  are  rounded  off. 

As  mentioned  in  Section  2.2,  except  for  the  Wiener  process  and 
the  M/M/1  queue,  the  value  of  <*>(s)  must  be  obtained  using  a search 
technique.  Our  final  results  use  the  quadratic  formula  to  determine 
co(s)  for  the  Wiener  process,  but  in  all  other  cases  (M/D/1,  M/E^/l 
including  M/M/1,  and  the  gamma  input  process)  we  use  the  standard 
Newton-Raphson  method.  This  search  technique  was  chosen  because  expres- 
sions for  the  derivatives  of  the  functions  were  available  and  subsequent 
computations  demonstrated  that  the  search  usually  converged  in  only  four 
or  five  iterations. 

The  exponent  function  $(•)  is  designated  FMD1,  FMEK,  or  FGAM 
in  the  FORTRAN  subroutines,  and  functions  F0MD1,  FOMEK,  and  FOGAM  are 
defined  as 

FO  •••  Ms)]  = F **•  [w(s)]  - s 
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Thus,  in  order  to  determine  'i)(s)  > 0 such  that  ^[^(s)  ] = s,  we 
search  for  the  value  of  u>  such  that  FO  •••  equals  zero.  Subroutine 
ZERO  performs  the  search,  computing  the  relative  accuracy  (change  in  co 
from  the  last  iteration) . When  the  absolute  value  of  the  relative 
accuracy  is  less  than  a specified  tolerance,  the  search  stops,  and  the 
current  value  of  u is  used  to  compute  the  Laplace  transform  using 
equation  (2.2.4). 

The  choice  of  tolerance  value  (FORTRAN  variable  TOLRNC)  affects 
the  accuracy  of  oo  determined  from  the  search,  and  the  accuracy  of  to 
affects  p(s)  and  the  estimate  of  E[W(t)].  Figure  3.5  shows  both 
the  approximations  of  E[W(t)]  based  upon  the  exact  quadratic  solution 
of  the  functional  equation  (2.2.4)  and  the  approximations  obtained  using 
the  Newton-Raphson  search  with  various  specified  tolerances.  We  observe 
that  the  two  results  agree  for  nine  or  ten  signficiant  digits  in  all 
cases  and  that  the  number  of  search  iterations  does  not  increase  greatly 
as  the  tolerance  is  decreased.  For  each  of  the  ten  epochs,  w(s)  is 
evaluated  thirty-four  times  (N  ■ 34);  thus,  the  total  number  of  itera- 
tions shown  in  Figure  3.5  applies  to  340  evaluations  of  w(s) . Since 
the  average  number  of  iterations  in  each  evaluation  only  increases  from 

4.2  to  5.3  as  the  tolerance  decreases  from  10  ^ to  10  we  have  set 

-20 

the  search  tolerance  at  10 

In  Chapter  2 we  discussed  the  stochastic  process  W(t)  in  an 

M/G/l  queue  and  its  expected  value  function.  Equation  (2.1.3)  is 

restated  here  in  a slightly  different  form  along  with  the  Laplace 

transform  of  E [W( t ) ] , equation  (2.2.4): 
z 
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E [W(t)]  = z - (l-p)t  + E [I(t)]  , 

z z 


i 


Vz’s)  - f 


-w(s)  z 
e 

s (s) 


The  first  two  terms  of  PyCz.s)  can  be  inverted  "by  inspection,"  i.e., 

by  referring  to  any  table  of  function-transform  pairs.  Only  the  last 

term  requires  numerical  inversion  using  LINV.  Stehfest  [1970,  p.  48] 

cautioned  the  prospective  user  of  LINV  as  follows:  "One  ought  to  be 

sure  a priori  that  the  unknown  function  F(t)  has  not  any  discontinuities, 

salient  points,  sharp  points,  sharp  peaks,  or  rapid  oscillations." 

Recall  that  our  unknown  function,  E [ I ( t ) ] in  M/G/l  queues,  is  a non- 

z 

decreasing  function  and  that  it  must  be  zero  between  t = 0 and  t = z. 

Although  it  is  not  a discontinuous  function,  its  slope  does  change 

abruptly  at  t * z,  and  we  might  expect  some  inaccuracies  near  that 

point.  The  top  section  of  Figure  3.6  shows  scaled  approximations  of 

[ I ( t ) ] evaluated  at  t * z for  M/M/1  queues  with  traffic  intensities 

between  .5  and  2 and  scaled  initial  server  loads  between  .2  and  4.  If 

the  inversion  technique  is  accurate,  then  all  entries  in  the  table  should 

be  zero.  We  observe  that  there  are  some  inaccuracies,  particularly  for 

queues  with  low  traffic  intensities.  For  example,  the  worst  case  is 

0 ■ .5  with  z'  * 1.0,  where  scaled  E [I(t)]  is  .007  instead  of  zero. 

z 

Since  the  inaccuracies  may  be  related  to  the  discontinuity  of 

the  derivative  of  E [ I ( t ) ] at  t - z,  our  approach  does  not  invert 

z 

the  third  term  of  Pw(z,s)  directly.  Instead,  we  construct  a new 
function. 
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3.6  CORRECTIONS  FOR  DISCONTINUOUS  FIRST  DERIVATIVE 
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where  g(t)  and  its  Laplace  transform  are  known,  and  where  g(t)  is 

chosen  so  that  both  H(t)  and  its  derivative  are  continuous.  Then 

LINV  is  used  to  invert  the  Laplace  transform  of  H(t),  and  approximate 

E [ I ( t ) ] is  obtained  by  subtracting  g(t)  from  approximate  H(t). 
z 

It  can  be  shown  that  the  slope  of  Ez[I(t)]  in  M/G/l  queues 

is  equal  to  the  probability  that  the  server  load  is  zero,  denoted 

P{W(t)  =0},  an  any  epoch  t.  Therefore,  the  slope  of  Ez[I(t)]  at 

t * z is  P{W(z)  = 0},  which  is  the  probability  that  no  arrivals  will 

occur  between  t = 0 and  t = z,  or  e . For  H(t)  to  have  a 

continuous  first  derivative  at  t = z,  that  derivative  must  be  zero. 

_\z 

Since  the  slope  of  E [ I ( t ) ] is  e at  that  point,  then  the  slope 

z 

of  g(t)  must  be  -e  ^Z.  Therefore,  we  define  g(t)  as  follows: 


for  0 < t < z 


g(t) 


-e  Xz(t-z)  , 


for  t > z 


Referring  to  any  table  of  transform-function  pairs,  the  Laplace  transform 


of  g(t)  is  -exp[-z(X+s) ]/s  . Thus  the  Laplace  transform  of  H(t), 


i.e.,  the  transform  of  E [ I ( t ) ] + g(t),  is 

z 


(3.3.1) 


» -oo(s)z  -z(\+s) 

PH(z,s)  - / e-S  H(t)  dt  « - * — 

U S 
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Applying  LINV  to  P (z,s)  we  obtain  an  approximation  of  H(t) . Using 

overscore  to  denote  approximations,  we  compute  the -estimate  of  E [ X ( t ) ) 

z 

as  follows: 


/ 0 for  0 < t < z 

Ez[I(t) ] = H(t)  - g(t)  = H(t)  + j 

' e ^Z(t-z)  for  t > z 

The  bottom  section  of  Figure  3.6  shows  scaled  E [I(t)]  when  the 

z 

correction  term,  g(t),  is  included  in  the  analysis.  The  results  are 

considerably  improved,  and  since  the  main  tables  will  express  E [W(t)] 

z 

in  hundredths  or  thousandths,  the  slight  inaccuracies  which  remain  will 

not  affect  our  tabulated  results. 

The  correction  term,  g(t),  is  employed  in  our  computations  for 

both  M/D/1  and  M/E^/l  queues.  Equation  (3.3.1)  is  incorporated  into 

subprogram  functions  PMD1  and  PMEK  of  the  main  programs.  The  gamma 

input  process  does  not  require  the  correction  term  because  both  E^ [ I ( t ) ] 

and  its  derivative  are  continuous,  i.e.,  the  slope  of  E [ I ( t > ] is  zero 

z 

at  t = z.  The  method  is  not  applied  to  the  Wiener  process  because 
process  W(t)  cannot  be  separated  into  four  components,  z + S(t)  - t + I(t) ; 
thus,  the  third  term  of  Pw<z,s)  cannot  be  interpreted  as  the  Laplace 
transform  of  expected  cumulative  idleness  in  this  case. 

So  far,  our  discussions  in  this  section  have  covered  some  of  the 
important  details  of  the  two  main  computer  programs  which  produce  the 

tables  of  scaled  E (W(t)l  is  Appendix  C.  Listings  of  these  two  main 
z 
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programs,  one  for  p = 1 and  one  for  p f 1,  are  included  in  Appendix  B. 

We  now  discuss  the  algorithm  for  M/E^/l  queues  with  p ^ 1,  and  we 

briefly  mention  the  differences  in  the  algorithms  for  the  other  three 

processes  and  for  p = 1.  Figure  3.7  describes  both  the  steps  of  the 

algorithm  and  the  names  of  the  subprograms  which  perform  those  steps. 

The  first  program  specifies  the  values  of  the  parameters  for 

subroutines  LINV  and  ZERO.  Values  of  p,  scaled  initial  server  load 

(z*),  and  sixteen  scaled  epochs  (t*)  are  read  from  punched  cards. 

The  program  performs  computations  for  three  pages  of  tables  on  each 

run.  The  three  tables  on  a page  have  the  same  value  of  p and  epochs 

t*,  but  each  table  has  a different  value  of  z*.  For  each  table  we 

compute  Ez[W(t)]  for  the  Wiener  process  at  all  sixteen  epochs,  then 

for  the  M/D/1  queue  at  the  same  sixteen  epochs,  followed  by  each  of 

the  seven  M/E^/l  queues  and  the  gamma  input  process. 

For  each  of  the  two  processes  the  first  step  is  conversion  of 

z*  and  t*  to  unsealed  values.  As  discussed  in  Section  2.3,  the 

2 2 

time  scale  conversion  factor  is  (l-p)  /cr  ; this  factor  is  denoted 

ALFSQR  in  the  subprograms  RZWNR,  RZMD1,  RZMEK,  and  RZGAM.  The  value 
2 

of  a for  the  unsealed  process  is  selected  so  that  the  formulas  and 

numerical  computations  are  subsequently  simplified.  The  conversion 

2 

factor  for  server  load  or  wait  is  |p.|/a  , denoted  WTFCTR  in  the  sub- 
programs. 

The  initial  trial  value  of  co(s)  in  the  Newton-Raphson  search 
for  the  solution  of  <H«(s) ] * s is  denoted  OSTART.  For  a specific 
z*  and  p < 1,  OSTART  is  arbitrarily  set  equal  to  1 for  the  first 
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FIGURE  3.7.  Flowchart  of  Algorithm  for  Tables 
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search  pertaining  to  the  first  epoch  t*.  In  preliminary  investigations 
this  first  search  usually  converged  in  seven  or  eight  iterations.  All 
subsequent  searches  related  to  the  specific  z*  and  p use  the  previous 
value  of  co(s)  as  OSTART,  usually  converging  in  four  or  five  iterations. 
For  p > 1 it  is  possible  that  the  Newton-Raphson  search  could  find  a 
solution  of  $[oo(s)]  = s such  that  co(s)  < 0.  To  ensure  that  the  search 
converges  to  co(s)  >0  in  this  case,  we  determine  OSTART  using  a two- 
stage  search  as  described  in  Figure  3.7. 

If  the  epoch  to  be  evaluated  is  less  than  or  equal  to  the  initial 
server  load  z,  then  the  expected  cumulative  idleness  is  zero.  For 
t > z we  use  subprogram  LINV  to  numerically  invert  the  Laplace  transform 
of  [ I ( t ) ] . On  any  given  run  of  the  computer  programs  (when  three 

pages  of  tables  are  prepared),  the  array  is  computed  only  when 

subprogram  LINV  is  first  called.  This  feature  and  others  are  explained 
in  a fully  documented  version  of  LINV  included  in  the  first  section  of 
Appendix  B,  "Computer  Program  for  Figure  3.2." 

The  last  section  of  subprogram  LINV  evaluates  the  Laplace  transform 
at  34  values.  Each  evaluation  requires  a Newton-Raphson  search  to 
determine  co(s) , except  for  the  case  of  the  Wiener  process  where  OMEGA 
is  determined  by  the  quadratic  formula.  The  number  of  iterations  in 
the  search  is  reduced  by  using  the  current  value  of  oo(s)  as  the  initial 
trial  value  for  the  next  search.  Subprogram  LINV  uses  array  V^  and 
the  34  values  of  the  Laplace  transform  to  calculate  H(t). 

Subprogram  RZMEK  first  combines  H(t)  with  the  correction  term 

to  obtain  E [ I ( t ) ) , and  then  it  computes  E [W(t)]  for  the  unsealed 
z z 


I 
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process.  In  the  final  step  E [W ( t ) ] is  converted  to  scaled  form  and 

z 

stored  in  the  three-dimensional  array  SCLWT.  When  the  computations 
for  three  tables  (each  with  ten  processes  evaluated  at  sixteen  epochs) 
are  completed,  subprogram  RZOUT  prints  five  copies  of  the  tables  with 
varying  margins.  Then  the  main  program  reads  punched  cards  specifying 
p,  z*,  and  t*  for  the  next  page  of  tables. 

All  computer  programs  listed  in  Appendix  B were  coded  in  IBM 
FORTRAN  IV  [IBM,  1971]  and  processed  on  an  IBM  370/168  with  the  high- 
speed-mult iply  feature.  The  average  computer  time  needed  to  evaluate 
scaled  [W ( t ) ] for  one  process  at  one  epoch,  i.e.,  a single  value 

in  the  tables  of  Appendix  C,  was  approximately  0.17  second. 
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CHAPTER  4 


NUMERICAL  RESULTS  FOR  EXPECTED  SERVER  LOAD 

4.1.  Using  the  Tables 

In  this  section  we  explain  how  the  tables  in  Appendix  C can  be 
used  to  obtain  time- dependent  mean  server  load  in  M/E^/l  queues.  We 
will  repeat  only  the  notation  and  formulas  from  the  previous  chapters 
that  are  required  for  using  the  tables;  the  details  of  the  underlying 
theory  and  numerical  method  will  not  be  discussed  here. 

The  system  being  studied  is  the  standard  M/E^/l  queue.  We  assume 
that  the  analyst  has  specified  the  value  of  the  mean  arrival  rate, 
denoted  and  has  determined  that  the  service  times  (random  variable 
S)  can  be  described  by  a gamma  or  Erlang  distribution.  If  an  actual 
operating  system  is  being  studied  and  data  are  available,  the  analyst 
can  use  the  methods  described  by  Hora  [1978]  or  Reinmuth  [1971]  to 
verify  that  the  input  process  is  Poisson.  Similarly,  a chi-square 
goodness-of-fit  test  can  be  used  to  compare  the  actual  distribution 
of  service  times  with  tabulated  values  of  the  incomplete  gamma  function 
[Pearson,  1922].  In  cases  where  the  service  times  of  an  actual  system 
do  not  exactly  fit  the  gamma  distribution,  the  analyst  still  may  wish 
to  use  this  theoretical  distribution  as  an  approximation. 

We  define  the  probability  density  function  for  this  two-parameter 
gamma  distribution  using  the  mean,  E(S),  and  the  Erlang  shape  parameter, 
k (not  restricted  to  integer  values),  as  follows: 
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f(t) 


[k/E(S)]k  .k-1  -kt/E(S) 
(k-'l)!  C 6 


t > 0,  k > 0 . 


The  standard  deviation  of  this  service  time  distribution  is  E(S)/\/k, 

so  the  appropriate  value  of  parameter  k can  be  determined  either 

2 

from  the  mean  and  second  moment,  E(S  ),  or  from  the  mean  and  variance, 
Var(S),  as  follows: 


(4.1.1) 


mnf. 

Var(S) 


Figure  4.1  shows  two  groups  of  gamma  service  time  distributions, 
where  the  special  case  of  the  exponential  distribution  (k  = 1)  is 
included  in  both  groups.  The  mean,  E(S),  of  each  distribution  in  the 
figure  is  equal  to  unity,  so  the  standard  deviation  is  1/Vk  in  each 
case.  The  top  chart  includes  density  functions  with  k * 2 and  k ■ 4, 
i.e.,  with  less  variance  than  the  exponential  distribution.  The  limiting 
case  as  k becomes  infinitely  large  is  the  degenerate  distribution 
with  constant  service  times,  corresponding  to  the  M/D/1  queue.  The 
bottom  chart  of  Figure  4.1  includes  density  functions  with  k = .5 
and  k ■ .2;  these  two  distributions  can  be  used  to  approximate  systems 
where  the  service  times  have  more  variability  than  the  exponential  case. 

The  stochastic  process  being  studied  is  server  load,  denoted  W(t). 

■ 

This  process  is  also  called  unfinished  work,  system  backlog,  or  virtual 

I ■ 

waiting  time,  the  latter  because  W(t)  represents  the  time  that  an 
arrival  at  epoch  t would  wait  before  beginning  service.  The  system 
may  begin  operation  at  epoch  t ■ 0 with  some  initial  backlog,  W(0)  * z. 
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Our  results.  In  scaled  form,  are  for  time-dependent  mean  server  load, 

denoted  E [W(t)].  This  operating  characteristic  is  appropriate  when 
z 

waiting  costs  associated  with  an  actual  system  are  directly  proportional 
to  the  waiting  times  experienced  by  the  customers. 

To  use  the  tables,  the  analyst  must  describe  the  system  by 
specifying  X,  E(S),  k,  and  z.  Some  preliminary  computations  are 
required  before  consulting  Appendix  C,  and  these  steps  are  outlined  in 
flowchart  form  in  Figure  4.2.  We  use  an  asterisk  superscript  to  denote 
the  scaled  values  which  appear  in  the  tables,  e.g.,  E .[W*(t*)].  The 
scaling  factors  are  WF,  which  is  applied  to  the  virtual  waiting  time 
(server  load)  values,  and  EF,  which  is  applied  to  the  epochs  (points  on 
the  time  scale).  Values  of  p,  z*,  and  k should  be  determined  before 
consulting  Figure  4.3. 

The  tables  in  Appendix  C contain  scaled  results  for  eight  M/G/l 
queues:  the  M/D/1  queue  and  M/E^/l  queues  with  k = .2,  .5,  1,  1.5,  2,  3, 
and  4.  If  the  value  of  the  shape  parameter  k determined  by  the  analyst 
is  not  exactly  equal  to  one  of  the  seven  tabulated  values,  either  the 
closest  tabulated  value  of  k should  be  chosen  or  simple  interpolation 
can  be  employed.  In  most  cases  the  interpolated  values  can  be  determined 
by  mental  calculation  directly  from  the  scaled  results  for  adjacent  k 
values  in  a table.  It  should  be  noted  that  for  Interpolation  purposes 
the  gamma  input  process  and  the  M/D/1  queue  correspond  to  M/E^/l  queues 
with  k * 0 and  k **  “,  respectively. 

In  Appendix  C there  is  a separate  table  for  each  combination  of 
p and  z*.  If  the  analyst  chooses  not  to  use  the  closest  (p,  z*) 
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Note:  EF  * Epoch  Factor  for  scaling 
WF  - Wait  Factor  for  scaling 

FIGURE  4.2.  Flowchart  Instructions  for  Using  the  Tables 
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combination  which  is  tabulated,  then  interpolation  between  two  or  more 
tables  will  be  required.  Figure  4.3  illustrates  the  cases  which  might 
be  encountered,  and  the  appropriate  methods  for  these  cases  are  dis- 
cussed below.  For  any  specific  (o,  z*)  combination  being  studied, 
the  objective  is  to  determine  pairs  of  t*  and  Ez^[W*(t*) 1 to  which 
the  scaling  factors,  EF  and  WF,  will  be  applied. 

Case  A.  No  interpolation  is  required.  There  is  a table  in 
Appendix  C for  this  specific  (p,  z*)  combination,  and  the  values 
for  scaled  epochs,  t*,  and  scaled  mean  server  load,  EzA[W*(t*)],  can 
be  read  directly  from  the  table. 

Case  B.  The  exact  value  of  z*  is  tabulated,  but  the  exact  value 
of  p falls  between  two  tabulated  values.  First  find  the  two  tables 
that  have  the  exact  z*  with  the  two  closest  values  of  p.  Then  list 
the  scaled  epochs  (t*)  that  appear  on  both  tables,  and  also  list  both 
values  of  ^z*[W*(t*)]  for  each  of  the  common  epochs.  For  each  epoch 
use  simple  interpolation  to  determine  Ez*[W*(t*)]  for  the  exact 
value  of  p. 

Case  C.  The  exact  value  of  P is  tabulated,  but  the  exact  value 
of  z*  falls  between  two  tabulated  values.  The  method  is  similar  to 
Case  B.  Find  the  two  tables  with  the  exact  P and  the  two  closest 
values  of  z*;  list  the  common  values  of  t*  with  the  two  values  of 
E .[W*(t*)];  and  use  simple  interpolation.  The  initial  behavior  of 

Z" 
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the  system  can  also  be  determined  as  follows: 


for  p < 1,  0 < t*  < (l-p)z* 

for  p = 1,  0 < t*  < z* 

for  p > 1,  0 < t*  < (p-l)z* 

Case  D.  The  exact  values  of  neither  p nor  z*  are  tahulated, 
but  the  exact  value  of  each  does  fall  between  two  tabulated  values.  For 

9 

one  of  the  closest  z*  values  which  are  tabulated,  perform  the  inter- 
polation between  the  two  closest  values  of  p using  the  procedure  of 
Case  B.  Use  the  same  procedure  for  the  other  closest  value  of  z*. 
Finally,  use  the  procedure  of  Case  C to  interpolate  between  those  two 
sets  of  data.  For  example,  if  we  are  interested  in  results  for  (p  = .67, 
z*  = .44),  we  interpolate  between  (p  = .6,  z*  = .4)  and  (p  = .7, 
z*  « .4)  to  obtain  (p  = .67,  z*  = .4) . We  also  interpolate  between 
(p  * .6,  z*  = .6)  and  (p  = .7,  z*  = .6)  to  obtain  (p  = .67,  z*  = .6). 

Finally,  we  interpolate  between  (p  = .67,  z*  = .4)  and  (p  = .67, 

z*  * .6)  to  obtain  results  for  (p  = .67,  z*  = .44). 

Case  E.  For  situations  where  p > 2 and  0 < z*  < 1,  the  Wiener 
process  provides  an  upper  bound  for  Ez*[W*(t*)]  in  any  queueing  system. 
Because  of  the  scaling  used  in  the  tables,  the  tabulated  values  of 
E JW*(t*)l  for  the  Wiener  process  depend  only  on  z*,  not  on  p. 

Z" 

Therefore,  any  table  with  p > 1.1  and  the  specified  value  of  z*  can 
be  used  to  determine  the  upper  bound.  Interpolation  can  be  used  if  the 


!z*  - t* 
z* 

z*  + t* 


exact  value  of  z*  falls  between  the  tabulated  values.  In  these  same 
situations  a lower  bound  for  Ez*[W*(t*)]  is  z*  + t*. 

Case  F.  For  situations  with  both  o > 1 and  z*  > 1,  a very 

accurate  approximation  of  Ez^[W*(t*)]  is  z*  + t*.  The  accuracy  is 
best  when  both  p and  z*  are  large,  but  even  for  p = 1.1  and 
z*  * 1 the  maximum  error  (M/D/1  queue,  t*  = 1.5)  is  only  1.4  percent. 

Case  G.  For  p = 1 and  z*  > 4,  the  initial  behavior  is 
EzJk[W*(t*)]  = z*  for  0 < t*  < z*.  For  t*  > z*,  the  results  for 
z*  * 4 provide  a loose  lower  bound. 

Case  H.  For  .5  < p < .9  and  z*  > 4,  the  initial  behavior  is 

Ez*[W*(t*)]  = z*  - t*,  for  0 < t*  < (l-p)z*.  For  t*  > (l-p)z*, 

the  results  for  z*  = 4 provide  a lower  bound. 

Case  I.  For  systems  with  p < .5  the  approach  to  equilibrium  is 
very  slow.  If  the  initial  server  load  z*  > 0,  the  initial  behavior 
is  the  same  as  Case  H.  For  z*  < 4 the  results  for  the  same  queue 
(i.e.,  the  same  value  of  k)  with  p = .5  provide  an  upper  bound. 

Other  Cases.  For  situations  with  .9  < p < 1.0  and  0 < z*  < 4, 
the  tabulated  results  for  the  Wiener  process  with  any  p < 1 and  the 
appropriate  z*  provide  a tight  upper  bound  for  any  queueing  system. 
Likewise,  for  situations  with  1.0  < p < 1.1  and  0 < z*  < 4,  the 
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tabulated  results  for  the  Wiener  process  with  any  P > 1 and  the 
appropriate  _z*  provide  a tight  upper  bound.  In  both  cases  the  Wiener 
approximation  is  most  accurate  for  values  of  p closest  to  unity. 

The  tables  in  Appendix  C can  be  used  along  with  Figures  4.2  and 
4.3  to  obtain  explicit  results  for  the  transient  behavior  of  mean 
server  load  in  a wide  range  of  M/G/l  queues.  In  many  applications 
the  analyst  will  perform  some  sensitivity  analysis  by  using  various 
estimates  of  X,  E(S),  k,  and  z and  comparing  the  results.  In  cases 
where  p < 1 the  analyst  may  be  interested  only  in  the  speed  with 
which  the  system  approaches  steady-state.  Expressing  the  Pollaczek- 
Khintchine  formula  with  our  notation,  the  mean  of  the  steady-state 
distribution  for  server  load  (often  called  expected  waiting  time  in 
the  queue,  excluding  service)  is 

E [W(->]  = WF/2  , z > 0,  p < 1 

z - 

4.2.  Several  Sample  Problems 

In  this  section  we  illustrate  the  methods  described  in  Section  4.1 
by  working  several  sample  problems.  The  first  three  problems  show  how 
changes  in  the  variance  of  service  times  affect  the  transient  behavior 
of  mean  server  load  and  the  steady-state  value,  assuming  the  other 
parameters  of  the  system  are  held  constant.  The  fourth  and  fifth  problems 
illustrate  using  interpolation  to  determine  results  for  queues  when  the 
exact  value  of  p or  z*  is  not  tabulated.  The  determination  of  waiting 
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costs  using  time-dependent  mean  server  load  is  discussed  and  applied  to 
one  of  the  sample  problems. 

The  following  problems  are  adapted  from  an  example  described  by 
Hillier  and  Lieberman  [1974,  p.  436].  A university  plans  to  lease  a 
small  batch-processing  computer  for  its  students  to  use.  It  is  expected 
that  the  students  will  submit  programs  to  be  run  every  3 minutes  on 
the  average  and  that  the  times  between  submission  of  programs  have  an 
exponential  distribution.  The  computer  under  consideration  could 
process  an  average  of  25  typical  student  programs  per  hour  if  it  were 
run  continuously.  We  will  examine  the  transient  behavior  of  this 
M/G/l  system  for  an  eight-hour  period,  assuming  that  the  system  begins 
operation  with  no  backlog.  Thus  far  we  have  specified  X,  E(S) , and  z. 
For  Problem  A we  assume  that  the  processing  times  for  the  students' 
programs  are  exponentially  distributed,  i.e.,  k = 1.  Referring  to 
Figure  4.2,  we  perform  the  initial  computations,  using  hours  as  the 
time  unit. 

Problem  A. 

20,  E(S)  = 1/25,  k = 1,  z = 0 

20- (1/25)  = .8 

20‘2* (1/25) 2 ,, 

1-.8 

.32/. 2 = 1.6 
0/.32  = 0 


X = 
P = 
WF  = 
EF  = 
z*  = 


Referring  to  Figure  4.3  with  p = .8  and  z*  = 0,  we  see  that  Case  A 
is  appropriate,  i.e.,  the  exact  values  of  p and  z*  appear  in  a 
table  in  Appendix  C.  Referring  to  the  appropriate  table,  we  note 
that  the  sixteen  scaled  epochs  t*  are  between  .02  and  4.  Eight 
points  should  be  sufficient  to  describe  the  queue's  behavior;  the 
selected  values  of  t*  and  ®z*[W*(t*)]  from  Appendix  C are  shown 
below  in  the  first  two  columns.  Referring  to  Figure  4.2,  we  compute 
t and  E [W(t)]  by  multiplying  the  scaled  values  times  the  scaling 
factors.  The  desired  results  are  shown  below  in  the  last  two  columns. 


Problem  A. 


t* 

E .[W*(t*)] 
z* 

t 

E (W(t)] 
z 

.1 

.164 

.16 

.052 

.2 

.230 

.32 

.074 

.4 

.306 

.64 

.098 

.6 

.351 

.96 

.112 

1. 

.405 

1.6 

.130 

2. 

.461 

3.2 

.148 

3. 

.482 

4.8 

.154 

4. 

.491 

6.4 

.157 

The  steady-state 

mean  server 

load  in 

Problem  A is  WF/2=  .16 

hour  (9.6  minutes),  or 

four  average 

service 

times.  After  approximately 

one  hour  of  operations 

(t  * .96),  mean  server  load  (.112  hour)  in  this 

M/M/1  queue  is  about  70%  of  steady-state;  after  three  hours  it  achieves 
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90%  of  the  steady-state  value.  For  queues  with  p < 1,  EzJk[W*(t*)] 

approaches  a steady-state  value  of  .5;  therefore,  the  percentage  of 

steady-state  achieved  at  any  epoch  can  be  determined  simply  by  doubling 

the  E . (W*(t*)]  value, 
z* 

For  Problem  B we  consider  a system  with  the  same  arrival  rate 
and  the  same  mean  service  time,  but  we  assume  that  the  service  times 
are  constant  (deterministic) . In  the  context  of  our  original  computer- 
leasing problem,  it  might  be  that  the  student  programs  are  actually 
uniform  or  that  this  specific  computer  has  very  little  variation  in 
processing  times  for  the  programs  being  run.  Referring  to  Figure  4.2, 
we  perform  the  computations  for  this  M/D/1  queue. 

Problem  B. 

X » 20,  E(S)  = 1/25,  k = °°,  z = 0 
p - 20.(1/25)  * .8 
WF  - . .16 

EF  = .16/. 2 = .8 
z*  = 0/.16  * 0 

Referring  to  Figure  4.3,  we  see  that  the  exact  values  of  p and  z* 
are  tabulated,  and  we  select  eight  points  from  the  appropriate  table 
in  Appendix  C as  shown  below.  The  tabulated  values  are  mutiplied  by 
EF  and  WF  to  obtain  t and  Ez[W(t)],  respectively. 


Problem  B. 


t* 

Ez*[W*(t*) J 

t 

Ez(W(t) ] 

.2 

.242 

.16 

.039 

.4 

.316 

.32 

.051 

.6 

.360 

.48 

.058 

1. 

.412 

.8 

.066 

1.5 

.446 

1.2 

.071 

2. 

.465 

1.6 

.074 

3. 

.484 

2.4 

.077 

4. 

.492 

3.2 

.079 

The  steady-state 

mean  server 

load  in 

Problem  B is  WF/2  = 

.08 

hour  (4.8  minutes),  or 

two  service 

times. 

We  observe  that  after 

a 

half-hour  of  operation 

(t  = .48)  mean  server  load  (.058  hour)  is 

about 

70%  of  steady-state;  after  one  and  a half  hours  (t  = 1.6),  It  has 
reached  90%  of  the  steady-state  value  (.074  relative  to  .08).  Compared 
with  the  M/M/1  queue  of  Problem  A,  this  M/D/1  queue  has  a lower  steady- 
state  mean  server  load  and  it  approaches  that  equilibrium  value  much 
faster. 

For  Problem  C we  consider  a system  with  the  same  prameters  as 
Problems  A and  B,  except  that  the  service  times  have  wide  variability. 
We  assume  that  the  mean  service  time  is  2.4  minutes  (1/25  hour),  as 
before,  but  that  the  standard  deviation  of  the  service  times  is  5.4 
minutes.  Referring  to  equation  (4.1.1),  the  appropriate  Erlang  shape 
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parameter  is  k = .2.  The  computations  from  Figure  4.2  and  the  selected 
points  from  Appendix  C are  shown  below. 


I 


Problem  C. 


X ss 

20,  E(S)  = 1/25, 

k * . 

2,  z * 0 

p = 

20- (1/25)  = .8 

WF  = 

20*6* (1/25) 2 _ 

96 

1-.8  = ‘ 

EF  = 

.96/. 2 * 4.8 

z*  = 

0/.96  = 0 

t* 

E ,*[W*(t*)] 

t 

Ez[W(t)] 

.1 

.156 

.48 

.150 

.2 

.223 

.96 

.214 

.4 

.299 

1.92 

.287 

.6 

.345 

2.88 

.331 

.8 

.377 

3.84 

.362 

1. 

.400 

4.8 

.384 

1.5 

.438 

7.2 

.420 

2. 

.459 

9.6 

.441 

The  steady-state  mean  server  load  in  Problem  C is  WF/2  * .48  hour 
(28.8  minutes),  or  twelve  average  service  times,  which  is  much  higher 
than  the  M/M/1  and  M/D/1  cases.  We  observe  that  about  70%  of  steady- 
state  is  achieved  after  three  hours  and  that  approximately  eight 
hours  of  operation  are  needed  to  reach  90 7.  of  the  steady-state  value. 
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Figure  4.4  shows  mean  server  load  curves  for  Problems  A,  B,  and 
C,  where  these  three  queues  have  the  same  values  of  X and  E(S).  We 


observe  that  the  queues  with  higher  variance  of  service  times,  i.e., 
with  lower  values  of  k,  have  a higher  steady-state  mean  server  load. 
Furthermore,  the  approach  to  equilibrium  is  slower  for  the  queues  with 
the  higher  steady-state  value.  We  also  note  from  Figure  4.4  that  the 
ordering  of  the  curves  is  the  opposite  of  the  ordering  of  the  scaled 
values  in  the  table  (p  = .8,  z*  « 0)  in  Appendix  C;  this  reordering 
is  explained  by  the  different  scaling  factors,  WF  and  EF,  for  the 
three  queues. 

Returning  to  our  computer-leasing  problem,  we  next  consider  a 
computer  that  can  process  an  average  of  30  programs  per  hour,  i.e., 

E(S)  ■ 1/30.  We  assume  that  the  processing  times  are  exponentially 
distributed,  i.e.,  k ■ 1,  and  that  the  arrival  pattern  is  the  same  as 
Problems  A,  B,  and  C.  For  Problem  D we  will  assume  that  students  can 
leave  programs  during  the  night  for  processing  when  the  facility  opens. 
It  is  estimated  that  "about  a dozen"  programs  will  arrive  during  a 
typical  night;  thus,  the  initial  server  load  (assuming  12  jobs  at  1/30 
hour  each)  is  .4  hour.  Referring  to  Figure  4.2,  we  perform  the  follow- 
ing calculations. 
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[W(t) ] hours 


FIGURE  4.4.  Graphs  for  Sample  Problems  A,  B,  and  C 


Problem  D. 


X 

P 

WF 


20,  E(S)  = 1/30,  k = 1,  z = .4 

20- (1/30)  = .67 

20-2- (1/30)2  _ 

1-.67  " * 


EF  - .133/. 333  = .4 
z*  = .4/. 133  * 3 


Referring  to  Figure  4.3  with  p = .67  and  z*  = 3,  we  see  that  Case  B 

is  appropriate.  The  two  relevant  tables  are  (p  = .6,  z*  = 3)  and 

(p  * .7,  z*  = 3),  and  we  list  the  values  of  t*  appearing  on  both 

tables.  Those  thirteen  values  and  the  corresponding  values  of 

E Jw*(t*) ] from  the  two  tables  are  shown  in  the  first  three  columns 
z* 

below.  For  each  scaled  epoch  we  interpolate  to  find  the  approximate 

value  that  is  two-thirds  (or  about  seven-tenths)  of  the  difference 

between  the  values  for  p * .6  and  p = .7.  The  interpolated  value 

is  shown  below  in  the  fourth  column.  Finally,  we  obtain  t and 

[W ( t ) ] by  multiplying  t*  and  the  interpolated  values  of 

E .[W*(t*)]  times  EF  and  WF,  respectively, 
z* 

The  steady-state  mean  server  load  in  Problem  D is  WF/2  * .067 
hour  (4  minutes),  or  two  mean  service  times.  We  observe  that  transient 
mean  server  load  is  within  10%  of  the  steady-state  value  after  approxi- 
mately two  and  a half  hours  of  operation.  Figure  4.5  shows  the  curve 
for  Problem  D and  the  curve  for  the  same  queue  with  no  initial  load. 

(The  calculations  for  the  latter  case  are  not  shown  here.)  The  system 
with  z = 0 reaches  90%  of  steady-state  after  only  .8  hour  of  operation. 
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[ W( t ) ] hours 


FIGURE  4.5.  Graph  for  Sample  Problem 


For  our  last  sample  problem  we  consider  a computer  that  can  process 
an  average  of  only  20  programs  per  hour,  and  we  continue  to  use  the 
same  arrival  pattern  that  was  specified  for  the  first  four  problems, 
i.e.,  an  average  of  20  arrivals  per  hour.  Since  the  arrival  rate  is 
equal  to  the  service  rate  in  this  system,  the  traffic  intensity  parameter 
equals  unity  and  we  have  an  unstable  queue.  For  Problem  E we  will 
assume  exponentially  distributed  service  times  and  an  initial  server 
load  equivalent  to  five  average  jobs.  The  initial  calculations, 
following  the  steps  in  Figure  4.2,  are  shown  below. 

Problem  E. 

X = 20,  E(S)  = 1/20,  k ■ 1,  z * .25 
P » 20* (1/20)  = 1 
WF  = EF  = 2/20  = .1 
z*  = .25/. 1 = 2.5 

Referring  to  Figure  4.3  with  p = 1 and  z*  * 2.5,  we  see  that  Case  C 
is  appropriate.  From  the  tables  for  p = 1 in  Appendix  C we  note  that 
the  scaled  epochs  are  identical  in  the  tables  for  z*  * 2 and  z*  = 3. 

We  select  eight  of  the  sixteen  epochs  and  record  the  tabulated  values 
below.  From  the  discussion  of  Case  C in  Section  4.1  we  also  know 
that  Ez*[W*(t*)]  = 2.5  for  0 < t*  < 2.5. 


0 to  2.5 


2.5 


0 to  .25 


.25 


3 

2.056 

3.000 

2.528 

.3 

.253 

4 

2.149 

3.016 

2.583 

.4 

.258 

6 

2.362 

3.108 

2.735 

.6 

.274 

8 

2.577 

3.238 

2.908 

.8 

.291 

10 

2.784 

3.383 

3.084 

1. 

.308 

20 

3.680 

4.115 

3.898 

2. 

.390 

40 

5.052 

5.364 

5.208 

4. 

.521 

80 

7.068 

7.289 

7.179 

8. 

.718 

If  the  M/M/1  system  described  in  Problem  E continues  operating 
indefinitely,  the  mean  server  load  will  increase  without  bound.  In 
many  actual  systems  the  arrival  rate  equals  or  exceeds  the  service 
rate  for  short  periods  of  time,  but  such  situations  may  be  acceptable 
if  the  cost  of  providing  faster  service  exceeds  the  cost  associated 
with  the  high  waiting  times. 

The  five  sample  problems  show  that  the  methods  described  in 
this  chapter  can  be  used  to  determine  time-dependent  mean  virtual 
waiting  time,  an  important  operating  characteristic  of  M/G/l  queues. 
In  actual  applications  the  analyst  may  wish  to  determine  the  cost 


associated  with  E [ W ( t) ] when  analyzing  or  designing  a system.  In 
z 

their  discussion  of  decision  models  for  queueing  systems,  Hillier  and 
Lieberman  [1974,  Ch.  10]  discuss  appropriate  cost  functions  for  various 
problems  and  apply  their  methods  to  steady-state  distributions  of  both 
the  number  of  customers  in  the  system  and  the  waiting  times  of  individual 
customers.  As  previously  mentioned,  our  results  for  expected  server 
load  are  appropriate  for  evaluating  costs  only  in  systems  where  the 
waiting  costs  are  directly  proportional  to  waiting  times.  We  next 
describe  a method  for  evaluating  the  waiting  cost  in  situations  where 
the  transient  behavior  of  the  system  is  important. 

Let  c be  the  cost  per  unit  of  waiting  time  experienced  by  an 
individual  customer.  The  cost  associated  with  a customer  arriving 
at  epoch  t is  c*W(t),  where  there  may  be  some  initial  server  load 
W(0)  = z for  process  W(-)-  For  any  specific  epoch  t,  W(t)  is  a 
random  variable;  in  our  situation  the  cost  function  is  linear,  and  the 
expected  cost  for  a customer  arriving  at  epoch  t is  c*Ez[W(t)].  In 
cases  where  the  cost  function  is  non-linear,  we  would  need  more  informa- 
tion about  the  distribution  of  W(t),  not  just  its  mean,  in  order  to 
determine  the  expected  cost  of  a customer  arriving  at  epoch  t.  The 
probability  that  a customer  will  arrive  during  interval  (t,  t+At) 

is  Xdt,  and  such  a customer  would  incur  expected  cost  c*E  [W(t)]. 

z 

Thus,  the  expected  waiting  cost  during  a period  from  t = 0 to  t = T 
is 

T T 

(4.2.1)  / c E [W(t)]  Xdt  - c X / E [W(t)]  dt  . 

0 2 0 2 
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Note  that  the  integral  on  the  right-hand-side  of  (4.21)  is  simply  the 
area  under  the  curve  for  transient  mean  server  load. 


i 


I 


For  a queueing  system  that  approaches  steady-state  quite  slowly, 

the  analyst  should  use  the  integral  expression,  not  the  equilibrium 

mean  waiting  time,  to  evaluate  expected  cost.  To  illustrate  the 

difference,  consider  again  the  system  described  in  Problem  C.  Assume 

the  university  has  determined  that  for  this  situation  the  appropriate 

cost  associated  with  making  a student  wait  for  a program  to  be  processed 

is  $10  per  hour.  If  we  approximate  the  mean  waiting  time  using  the 

steady-state  value  of  .48  hours,  then  the  expected  waiting  cost  for  an 

individual  student  arriving  with  a computer  program  to  be  run  is  $4.80. 

The  mean  arrival  rate  is  20  per  hour,  so  the  expected  waiting  cost  is 

$96  for  each  hour  the  system  is  in  operation.  Thus,  during  an  eight- 

hour  day  the  total  expected  waiting  cost  is  $768.  This  total  cost 

figure  is  based  on  the  equilibrium  mean  waiting  time,  but  from  Figure  4.4 

we  can  see  that  E [W(t) ] for  Problem  C is  substantially  less  than 
z 

the  equilibrium  value. 

For  a more  accurate  determination  of  total  waiting  cost  during 

the  eight-hour  period,  we  can  evaluate  the  integral  on  the  right-hand- 

side  of  equation  (4.2.1)  graphically  by  plotting  the  curve  from  Figure  4.4 

on  a fine  grid  and  counting  the  squares.  Using  this  method  the  value 

2 

of  the  area  under  the  curve  is  approximately  2.66  hours  . (The  correspond- 

2 

ing  value  using  the  equilibrium  mean  server  load  is  .48-8  = 3.84  hours  .) 
Multiplying  the  area  under  the  curve  by  c-X,  we  find  that  the  correct 
total  expected  waiting  cost  for  the  eight-hour  period  is  $10‘20*2.66 
=■  $532. 
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In  this  particular  problem  the  use  of  the  equilibrium  value 

I I ' 

produces  a total  cost  figure  approximately  44%  greater  than  the  correct 
amount.  The  example  illustrates  the  importance  of  using  time-dependent 

l 

results  when  the  approach  to  equilibrium  is  slow.  The  graphical 
technique  can  also  be  used  to  determine  expected  waiting  cost  for 
unstable  queueing  systems  (p  > 1)*  where  all  behavior  is  transient 

| 

and  no  equilibrium  is  ever  achieved. 
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CHAPTER  5 
CONCLUSIONS 


The  objective  of  this  research  was  to  provide  tables  for  time- 
dependent  mean  server  load  in  M/G/l  queueing  systems.  Chapter  2 pre- 
sented the  theoretical  framework  for  the  server  load  process,  the 
Laplace  transform  of  Ez[W(t)],  and  a scaling  procedure  that  allowed 
us  to  reduce  the  number  of  parameters  required  to  describe  a specific 
system  from  four  [X,  E(S),  z,  k]  to  three  [p,  z*,  k] . In  Chapter  3 
we  considered  a technique  for  inverting  the  Laplace  transform,  and 
we  conducted  extensive  checks  and  comparisons  to  ensure  the  accuracy 
of  our  numerical  results.  The  sample  problems  in  Chapter  4 illustrated 
that  a practitioner  can  easily  use  our  tabulated  results  in  Appendix  C 
by  following  simple  step-by-step  procedures.  We  now  discuss  some 
additional  results  of  this  research,  including  a study  of  the  error 
associated  with  the  Brownian  approximation,  and  we  offer  some  suggestions 
for  future  research. 

The  scaling  procedure  employed  in  this  research  produces  curves 
for  Ez^[W*(t*)]  that  are  monotonically  ordered  by  each  of  the  three 
parameters  (p,  z*,  and  k) , and  these  consistent  orderings  facilitated 
the  interpolation  procedures  described  in  Chapter  4.  Referring  to  any 
single  table  in  Appendix  C,  i.e.,  fixing  the  values  of  p and  z*,  we 
observe  that  the  curves  for  Ez^[W*(t*)]  increase  monotonically  as  k, 
the  shape  parameter  of  the  Erlang  service  time  distribution,  increases. 
(This  ordering  of  the  transient  curves  by  the  shape  parameter  has  not 
been  proven  analytically  and  is  a possible  topic  for  future  theoretical 
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research.)  As  we  would  expect,  the  upper  bound  (corresponding  to 

k = “>)  is  the  curve  for  the  M/D/1  queue.  But  our  numerical  results 

also  indicate  that  a lower  bound  can  be  identified:  the  gamma  input 

process  (corresponding  to  k = 0) . We  can  think  of  these  upper  and 

lower  bounds  as  defining  an  "envelope",  and  we  might  hypothesize  that 

this  envelope  contains  the  curves,  in  scaled  form,  for  all  M/G/l  queueing 

systems,  not  just  for  the  M/E^/l  queues  that  are  studied  here.  Figure  5.1 

shows  the  upper  and  lower  bounds  for  the  (p  = .5,  z*  = 0)  case.  We 

observe  that  the  envelope  is  relatively  narrow;  specifically,  its 

maximum  width  is  11%  of  the  steady-state  value  in  this  case.  For 

z*  = 0 and  p = .6,  .7,  .8,  and  .9,  the  maximum  width  is  9%,  7%,  5%, 

and  3%,  respectively.  We  could  use  these  bounds  to  obtain  approximate 

2 

results  for  any  M/G/l  queueing  system  for  which  X,  E(S) , E(S  ) , and 
z are  known,  without  specifying  the  exact  shape  cf  the  service  time 
distribution. 

Another  monotonic  ordering  is  observed  by  fixing  the  values  of 
k and  z*  and  varying  p.  For  example,  if  we  examine  the  curves  for 
the  M/D/1  queue  (k  = °°)  with  z*  = 0,  we  observe  that  for  p < 1 the 
scaled  curves  increase  monotonically  as  p increases.  (In  Appendix  C 
this  comparison  requires  examining  tables  that  are  each  three  pages 
apart.)  The  M/D/1  curves  seem  to  be  approaching  the  curve  for  the 
Wiener  process  as  p increases.  This  observation  has  not  been  proven 
analytically,  but  the  heavy  traffic  (or  Brownian  motion)  approximation 
for  the  equilibrium  distribution  of  M/G/l  server  load  is  derived  by 
examining,  the  scaled  form,  the  limit  of  a sequence  of  queues  as  p 
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research.)  As  we  would  expect,  the  upper  bound  (corresponding  to 
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lower  bounds  as  defining  an  "envelope”,  and  we  might  hypothesize  that 
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for  the  equilibrium  distribution  of  M/G/l  server  load  is  derived  by 
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FIGURE  5.2.  Percentage  Error  of  the  Wiener  Approximation,  z*  * 0 
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relationship  using  a quadratic  function.  In  such  cases  knowledge  of 
both  the  first  and  second  moments  of  the  server  load  distribution  is 
required.  Prabhu  [1965]  derived  the  double  Laplace  transform  for  the 
time-dependent  server  load  distribution  in  M/G/l  queues.  If  we 
evaluate  the  second  derivative  of  that  expression  at  zero,  we  obtain 
the  Laplace  transform  for  the  time-dependent  second  moment  of  the 
distribution.  The  LINV  algorithm  can  then  obtain  the  second  moment 
at  specified  epochs  for  use  in  a quadratic  cost  function.  There  may 
be  other  queueing  applications  for  this  relatively  efficient  and 
accurate  numerical  technique  for  Laplace  transform  inversion. 
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It 


approaches  unity.  Our  empirical  observations  suggest  that  for  both 
the  transient  and  equilibrium  cases  it  might  be  possible  to  obtain 
analytic  proofs  for  the  ordering,  not  just  for  the  limit. 

The  Wiener  process  was  included  in  our  research  because  it  is 

often  proposed  as  an  approximation  of  server  load  in  M/G/l  queues. 

If  our  detailed  results  for  M/E^/l  queues  are  not  available,  an 

analyst  can  use  the  scaling  factors  and  a brief  tabulation  of  the  scaled 

Wiener  process  to  obtain  such  an  approximation.  From  Figure  5.1 

we  can  see  that  the  method  will  be  most  accurate  for  M/D/1  queues  or 

M/E  /I  queues  with  large  values  of  k.  Conversely,  the  error  will  be 
k 

greatest  when  the  method  ia  used  for  an  M/E^/l  queue  with  k near  zero. 
The  limit  of  this  "worst  case"  is  the  gamma  input  process,  and  our 
numerical  results  in  Appendix  C allow  us  to  make  a detailed  study  of 
the  possible  error  associated  with  the  Brownian  approximation.  We 
express  the  difference  between  the  values  for  the  Wiener  approximation 
and  the  gamma  input  process  as  a percentage  of  the  value  for  the  gamma 
input  process,  and  Figure  5.2  tabulates  this  error  for  all  values  of 
P at  selected  scaled  epochs.  We  observe  that  this  "worst  case"  error 
decreases  monotonically  as  p approaches  unity  or  as  t*  increases. 

This  research  used  numerical  inversion  of  the  Laplace  transform 
to  determine  the  first  moment  of  the  server  load  distribution  at 
specified  epochs.  In  Chapter  4 we  used  the  transient  curve  to  deter- 
mine waiting  cost  for  cases  where  the  cost  is  directly  proportional 
to  the  customer  waiting  time.  In  many  applications,  the  cost  function 
may  not  be  linear,  and  the  practitioner  may  approximate  the  cost 
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APPENDIX  A 


SUMMARY  OF  NOTATION 

This  list  includes  brief  definitions  of  the  notation  used  most 
frequently  in  this  dissertation. 

A(t)  stationary  Poisson  process,  the  M/G/l  arrival  process 

c cost  per  unit  of  waiting  time 

E(*)  expectation,  first  moment,  of  a random  variable 

E [•]  expected  value,  conditional  on  initial  load  z 

z 

EF  scaling  factor  for  epochs 

Fa  an  approximation  based  on  Stehfest's  algorithm 

F(*)  cumulative  probability  distribution  for  S 

F*(s)  Laplace-Stieltjes  transform  of  F(-) 

f(')  probability  density  function 

fn(t;  a)  observational  density  function  with  parameters  a and  n 

G(t)  gamma  process 

g(t)  correction  term  function 

H(t)  Ez [I(t) ] + g(t) 

I(t)  cumulative  idleness  process 

k Erlang  shape  parameter 

N number  of  values  of  the  Laplace  transform  used  by  Stehfest's 

algorithm  to  obtain  an  approximation 

P(t)  any  function  of  interest 

P an  approximation  of  P(t)  based  on  n values  of  the  Laplace 

transform  using  Gaver's  method 
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• 

PH(z*s) 

" " ■ n 

Laplace  transform  of  H(t),  conditional  on  z 

Pw(z,s) 

Laplace  transforn  of  W(t),  conditional  on  z 

P(s) 

any  Laplace  transform  function 

j s’  si 

service  time  random  variable 

S(t) 

compound  Poisson  process,  an  input  process 

, T 

random  variable  for  an  observational  epoch 

Var[- ] 

variance  of  a random  variable  or  process 

W 

random  variable  for  equilibrium  server  load 

WF 

scaling  factor  for  waiting  time 

W(t) 

server  load,  or  virtual  waiting  time,  process 

X(t) 

net  input  process,  a Levy  process 

z 

deterministic  initial  server  load,  W(0) 

“i 

error  components  of  an  asymptotic  expansion 

a,  p 

shape  and  scale  parameters  for  gamma  density  function 

r(-) 

gamma  function 

X 

mean  arrival  rate,  parameter  for  A(t) 

-E[X(t)]/t 

C(t) 

standard  Wiener  process 

P 

traffic  intensity  parameter 

2 

CT 

Var[X(t) ]/t 

4>(-) 

exponent  function 

w(s) 

solution  of  functional  equation  4>[co(s)]  = s 

” 

overscore,  or  "bar",  denotes  approximation 

* 

asterisk  superscript,  denotes  a scaled  value 
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APPENDIX  B 


LISTINGS  OF  THE  COMPUTER  PROGRAMS 

1.  Figure  3.2,  including  the  documented  version  of  LINV 

2.  Figure  3.4 

3.  Figure  3.5 

4.  Figure  3.6 

5.  Appendix  Tables,  p 4 1 

6.  Appendix  Tables,  p = 1 
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CuflrOT£fc  £aOG)Uft  fOE  *J.UU8t  3.2 


lHkA.li.il  EE  Ai.*6(A-ri,0-2)  , IEaIGaE  (1-6) 
Oufc1ic6/i.A£LAC/ErOC8  (10)  , A?£  Eu  X (1 C ,0 ) 
olfatfcbiUA  EXAC1  (10,6)  ,51'EHE  (10,6)  , V(50) 
x.AX£a6A2  £1  ,*2  ,?3  ,£4  ,*5,£0 
ti-0 
6 = 34 

A.A1A  Sl£tt£  /.365j5,.  39912,. 32655,  .c6278, 

1.2  51  74,  .22  569,  .21322,.  1993b,.  16b  14,.  17796, 

1.  1to_>66,  1.  3^543,4.47354, 

110. b034 2, 2C. 70 64 5, 35. 7 6632, 36 .62 533, 64. o2735, 
1120.76473,165.66749, 

1 . 6 6775,  .91 001,  .O3826,.30  9oo, 

1-.  0 21 16,-.  31  62  7,  - .37254,  -.70609,-.  9 1049,-.  56949, 

1-.  37  /62,-1.  27084,-1.67544, 

1-1 .90352,-2.16727 ,-2. 36670 ,- 2. 52270, -2 .o574C , 

1-2. 77390,-2  .6o 091 , 

1.  36756,  .13357,  .05C43,  .01845, 

1.CGo40, .00 195, .00036,- .OOOOo ,- .00047 ,-. uOO 20 , 
1-.cb333,-.32531,1.02575,2.39533, 

12.  76644,  1.21092, -3. 32956, -11.  a 29  53, -25  . 2 63  >3, -44, 663 11/ 
Lo  30  1=1, 1C 
££  UCH(i)=i 
1-1 

EXACX  (1,1)  = 1/J>5wa1  {1*3. 1415926335697532364626433) 

EXAcX  (1,2)  =X*l*T/6 

LX AC'l  (1,3)  =OiXli  (Oowax  (2*1)  ) 

EX  ACl  (1,4)  =-D2GG  ('x)-O.  57 7 2l36o 45 01  53 366061 
ilkCl  (1,5)  =D£XE(-1) 

EXAC'i.  (x,6)  = 1-  3*1  ♦ 3*1  *1/2-1*  1* l/o 
cAll  i,mk'  (t  1,5  ,1  , A££j\02  (1,1)  , V.ri) 

CALL  01 NV  (£>2,6,1,  A*>£aOX  (1,2)  , V,d) 

CAjliA.  2IaV  (£3,6,1,A££aGXII,3)  , V,tt) 

CA22  216V  (14,6  ,1  , a££EGX  (1,4)  , V,fc) 

CALL  4.I6V  (£5,6,1,  AtPEuX  (1,5)  , V , d) 

CALL  4.1 6V  ( £6,  6,1,  A£f*0X(2,6)  ,V,d) 

30  C0M1I6UE 
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maIIE  (6,4) 

4 60EMA1  (161,////) 

6AEG16  = 4 ♦ 6C0PX 
00  6 xdfiG6  = 5,dA&GlM 
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6 fCfctiAl  (1H  ) 

6 U>61160£ 

9EITt  (6  ,40) 

40  E0EHAl(l6  ,332,'AIGUEE  3.2  16VEE5IG6  0£  1£61', 

1'  FU  60*1063  ') 

ME  11  £ (6  ,42) 

42  £0aqA1‘  (1  h ,352,39  (*-•)  ) 

■ EITE  (6,43) 

43  AOaEAI  (1i»0,33X,2  ('APPECXIdAlE  E(l)  U5I6G*  ,21X)  ) 


v-ufl.  OTEP.  PaUGfaAh  ton  tibUdc  3.2 


IfiPi.Icil  6E  Ai.*blA-il,0-Z)  , lAX  iGoi  (I-N) 

CO  6dc6/i.A  PLAC/EPC  Cd  ( 10 ) , APiduX  (10  ,0  ) 

LitibbblOb  EXACT  (10,6)  ,5lEd)r  (10,6)  , 7(50) 

oXTEnAAL  PI  ,t2,P3  ,P4  ,*5,P0 

ft=0 

6=34 

bA'iA  S1EH*  /.3663b,.  35)91  A,.  3266b,  .*8278, 
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1.  16366,  1. 34.5*0,4. 47  354  , 

1 10 . to 034  2, 20 .70b45,3b. 76632, 36 . 62533 , b4.  o27  3b , 

IIaO. 76473, 165. 667 49, 

1. 9 6775,  .91001,  .b562b, .30 boo, 

1-.0  2115,-.  31  So 7,  - .37254, -.7060  9,-.  91049,.-.  9654  9, 

1-. 57762,-1.  27084,-1.67544, 

1-1 .9o392,-2.16777 ,-2. 36670 , - *. 57270, -2 .o5740 , 
1-2.77350,-2  .bo051  , 

1.  36756, . 13b57, .05 C43, .01845 , 

1.C0o40, .00195, *000  36,- .OuOOo ,- .00047  ,-.o00  20  , 

1 -. 66333,-. 57331,  1.0 757 5, 2*  39533, 

12.  76644,  1.7 1092, -3. 32956, -11.  6 29  53,-25.263  >5, -44. efabll/ 
Lu  30  1 = 1,  1C 
EPUC  ki  (1)  =1 

4 = 1 

EXAC1  (1,1)  =1/DSwa'i  (1*3.  14 159 26 33 5697 ‘5 32364c 2 64 3 3) 

EXACT  (1,2)  =X*T*T/6 

EXAC1  (1,3)  = 1)316  (Oo  wax  (2*T)  ) 

EX  AC  x (1,4)  =-DEOG  (X) -0.577  2l56o  490 153  3 66061 
EXACi  (1,5)  =DEXP(-1) 

EXACI (x,0) = 1- 3*1 ♦ 3*X *1/2 -1*1* 1/6 
L&1.4  1.16  V (1  1,6  ,'X  , APPttOX  (I,  1)  , V,rt) 

CALL  LI SV  ((>2,6,1,  At-PaO X ( 1,  7)  ,1/  ,il) 

CAm.  Ll6  V (P  3,  N ,X  , APPaOX  (1,3)  , V,M) 
cALl  LI6  V (P4 ,6  ,1‘  , aPPEGX  (1,  4)  , V,6) 

CALL  *I6V (P5,6 ,1, APPajX  (1,5)  ,7,d) 

CALL  J.I6V  (96,6,1,  ApPdOX  (j.,6)  , Y,d) 

30  CONI'  I6U£ 

Do  60  6 COP  X = 1 ,5 
■ aIIE  (6,4) 

4 FGBHA1  (1H1,////) 

6A6016  = 4 ♦ 6CCPX 
00  6 4.6606  = 5 ,d AfiGIN 

6 All  b (6,0) 

6 fCfcBAl  (1H  ) 

6 Co  61 16  0 E 
6fiITE(0  ,40) 

40  PG6ftAl'(l6  , 3sX, ' * iGUdE  3.2  16VE66I06  OP  1E51', 

1'  FU  60x1063') 
mb  111  (6,47) 
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MfclTL  (6,4b) 
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12(2X,14('-'))  ,2X,9('-')) 

Mb IX e (o,50) 

FOBtUl  <1H0,2©X,'F  (1)  = 1/Sw*x  (PI*T)  • ,24X, 

1 • F (1 ) = -C-LE(X)') 

MBIT x,  (6,52) 

FOBeUX  (lii  , Ibl  ,2  (41  (•  - ')  ,2X)  ) 

Du  So  1=1 , 10  i 

bBri‘E(0,S4)  EA-GCE  (i)  , (EXACX  (I  , J)  , A A FBCX  (I  , 0 ) , 

1&X £Di  (l,o),  J*1»4/3) 

iCRflAHIfa  ,10i,F4.  1«  2 (2 (2X,F  14*10)  ,2X,f9.5)) 

CObllbUE 

mEIX  E (6 ,60) 

FOKhAX  (IhO,  26X,'F  (X)  = (X*1*X) /6 • ,27 X , 
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MBIXE  (6, 52) 
i/0  b6  1=1,10 

ME11&  (b  ,54)  EPOCa(I),  (EXACX  (I,  J)  , AFFEOX  (1,  J)  , 
ISXEaF  (i,  J)  , 0=2,5  ,3) 

CObXIbUE 
ME  HE  (o,7G) 

lOeoAX  (1d0,26X,'  £ ll)  = SIE  (bwEXU*T)  ) ' ,x2X, 
l'F(X)  = 1-3bl*3*I«!l/2-X*T*X/S  ') 

MB  lit  (o  ,52) 

LG  7b  1=1,10 

NbIX E (6, 54)  EFOCrt  (i)  , (EXACX  (1,0)  ,AFIBCX  (l,u)  , 
1SXEEF  (i,J)  , 0 = 3,  6,3) 
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CObXxbUE 
oX  OF 
EEL 
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CCABOD/GliX  FUx/wU  AL  (10)  , EXP  MX  (10,8)  ,IXbloX<8) 
xlMEbSIOb  E (25) , A (25)  , E<25),  V(50) 


C. . ABGLAcbl 
C 
c 
C 
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FUbCXlub  r(S)  * LaPaACE  1 BASSE  OBE  OF  F(I) 
EC.  OF  VALUES  OF  S USED  TO  APFuUaIHATE  f (X) 
VALUE  AX  iaiCti  F (X)  IS  XO  EE  AFPaOXlflAI  EU 
ABFHOXLBAX £ f (I)  bASEO  Ob  E VALUES  OF  ?(S) 
ABBA*  OF  COEFFICIEbXS  'US £0  abfEAxEOLY) 
t Ob  SAL  FAbAM  El'EE  EXrLAIMiiD  BELO«. 


«kIl'E  (6,44) 

FORMAT  (Iti  ,352,2 (17HST£UF£ SI'S  a£TaOD,2o2 )) 

MRIIE46  ,4b) 

4 OEM  AX'  (Id  , 322,2  (2b(*-*)  , 1dX>) 

Mb  HE  (6  #46)  R , M 

FubMAX(1a  ,112,*1 *,6X,2(*EXACX  F(X)  • , 92, ' R = ' ,I2,*2 , 
1*R=1C*,7  X)) 
wfcXTL  (6,4b) 

40baAI(lH  ,1C2,' *,2(22,  14  (*-*))  , 22, »(•-'.)  , 

12(22, 14(*-*))  ,22,94*-*)) 

Mb  HE  (o,bQ) 

40EMA1  (lH0,2oI,'f  (1)  * 1/Swaa  (PI*X)  • ,242, 

1 * 4'  (X1 ) *-C-Lb(l)*) 

MEL  IX i>  (6,92) 

FOUrtalilU  ,162,2  (41  (•-•)  ,22)) 

Do  So  1=1,10  ( 

AblXE  (0, 64)  EPOCh  (i)  , (EXAC1  (1,J)  ,A4ibC2  (I,J), 
loXEai  (X,o)  , J = 1, 4 , 3) 

FORMAT  (Ifa  ,102,Fb.  1,  2 (2(22, 4 14.10)  ,2X,F9.S)  ) 

CORXlRUE 

«EIX  £ (6,60) 

FORMAT  4 160,262, *4  (X)  = li*t*T) /6  • ,272, 

1 *4  (1)  = £24  (-!’)  ') 

MBIZE  (6,62) 

oO  6 6 1 = 1 ,10 

Mblli.(6,b4)  EPuCb  (2)  , (EXACT  (I, J)  , APPEGX  (1 , J)  , 
ISlEaF  (I,J)  , J=2,5  ,3) 
corxirue 
MR  IX  £ (o,70) 

iOonAX  (1d0,262,*  4 (1)  = SIR  (Swbl  ( 2**)  ) * ,X22, 

1*4(1)  = 1-3*l*3*i«ll/2-Z*T*T/6  ') 

WRITE  (o  ,62) 

00  7o  1=1,10 

maITE  (6,64)  EPOCh  (I)  , (EXACT  (1,0)  , AF4B02  (l,u), 

1 jT E lit  (a, J)  , 0 = 3,  6,3) 
to  MX  III  Ut 
CORXIRUE 
oIOF 
Eli  L 


oUBaCLxIAE  LIRV(4,R,1, FA, V,b) 

IMPLICIT  aEAL*  b (A-h,  A,  L,G-L)  , laTEGEb  (1,J,M,R) 
COMMCA/ IRPUI/  EPOC  fc  (10)  , LAMBDA,  fcAPSVC,RHO,UT0L(8) 
COaaOM/oUiPUx/wUA£(1Q)  , £2P MX  (10,8)  ,ITRluI(8) 
oIMERSIGR  £ (26) , 4(25),  6(26),  V pO) 


(EXACT  41, J)  , APPROX (1,J)  , 


t. , A&UUMERl 
t 


4 UR  Clio  R r (S)  * LAPi-ACE  1RAASFORE  Ok  F(I) 
RC.  04  VALUES  OF  S USED  TO  APBuGaIHATE  F (X) 
VALUE  AT  Malta  F (Z)  IS  TO  E£  APPaGXIMAI ED 
APPROX aMAX £ 4(1)  BASED  OR  R VALUES  OF  P(S) 
AREA*  04  COEFFICIENTS  'USED  oopEAlEDLY) 
k OR  UAL  PARAMETER  EIPLAIlUD  BELOm. 


■1 

IF  (M.LwM)  00  TO  bO 

1 

w 

C • •OtM 

FifcbL  CALL  OF  L1MV  (B  MOT  E^oA*  lu  M ) AaaAl  V(L>  filial 

I 

1 

C bt 

c 

EVALUATED. 

C..F**LT  CONFUTE  £ (I)  * FACTQaI AL  Of  21  oVE£  iuO bt  OF  1 kbit 

0 1-1 

, AML  A (I)  * FACTORIAL  OF  1-1. 

c 

MH  * H/* 

Ml)  * 2 
£ (2)  = 12 

F(1)  » 1 

*(2)  * 1 

1 

c 

DO  10  1 * 2,  M H 

£ (It  1)  * ( 2.D0*(2*I*1)*E  (I)  ) / I 

F(IO)  « F(I)*I 

I 

10 

C0MT1MUE 

c 

1' 

C..MEA1  Co&kll  JiL  ti  (*J)  * ALL  TLttt S IMOt?  LmDSMI  of  laLtl  1 la 

C ILL  V U)  LUHtlATlOM. 

C 

EMA  * Mti 

1 

C 

h ( 1)  * 2/f  (ME) 

1 

LO  20  J = 2,  Mti 

H (Jt  * t (*)•*»’ MB)  * £ (0)  ) /F  (Mh-J  ♦ 1) 

20 

cOMXIMUE 

c 

C..FIMALLX  EVALUATE  AttAAX  V (I)  . 

1 

C 

ILIUM  * 2*  (MH-2*  (a h/2)  ) - 1 

C 

LO  40  I s 1f  1 

V(I)  * 0 

JFIfiLX  * (1*1) /2 

OLALI  * hlM0(I,Mh) 

C 

LO  30  0 - JFLaLI,  JLASl 

V (I)  * V(I)  ♦ H(0)/(F(1-J  + 1)*F  l2*J-I*1)) 

30 

CO  MI  I MO  £ 

C 

VCD  * ILIUM* V (1) 

ILIUM  = -1SIUM 

40 

CCMTIMUE 

C 

fl  * M 

t 

C.  .SUbSEwUfcW  CAa.1l  (K  ti^OAL  TO  M)  Wii.L  USE  AaaAK  V (1)  FBU& 

C COB&OM  LTGAAOE  AMD  a ILL  LEAMCtt  '10  1U1S  SECIIOM  OF  LI  MV. 

0 

Eli  . 

50 



FA  * 0 

A * .6»3l471e0bb994b3094172J21t14bOD0  / 1 

88 

, 



£b  60  i«  1,  k 
kk  « F A ♦ V (I 
00  C0NX1NUE 


i A * A*£  A 


fcElUEN 

xNC 


DuUtoLE  fM£C IblON  EUNCIIOII  £>  1 <S) 
IWkjuICIl  REAx*8  (A-b,0-Z)  , INIEGeR(I-M) 
«1*1/D£M&1(S) 
kEXUE* 

£N  b 


DCUfliE  EBECISIO*  IfUNCTIO*  PE  (a) 
lapLiCiX  *£AL*e(A-Uro-Z)  , I*1EG£E  (I-N) 
t2=1/(b*J»*.»*i>) 

££'lb  £N 
tkb 


OOUbEE  P&ECIalON  JUMCXION  P3(a) 

Ifttx^Cxl  £EAL*8(A-ti,o-Z)  , INIEGEN  (I-li) 
P3=law*X(3.  141i*2  £535697  93  2384o2o433/  (2*S*a*S)  j 
1*o£iP(-1/(2*S)  ) 
aETUEN 
kJlb 


O OUtoxE  fcBECiaiGN  £ UN  Cl  10 II  £4(8) 
IfiPLlClI  Ei,AL*8(A-H,G-Z)  , Ifcl&GEK  (I-N) 
f4  =DxwG  (S)  /i> 

NEXUBN 

END 


OOUbxE  Pa*C1£ION  £ UNCI  ION  Pb(N) 

lEPIICIT  R£Ai.*8(A-h,u-Z)  , iNXEGEk  (I-N) 

rb- 1/  (S  + 1) 

kEiUNk 

END 


00 Udx£  PRECISION  HJNCTIcM  p 0 (a) 
IHPLICIX  k£Ai.*6(A-H#0-Z)  , IfcxtOEA(I-N) 
P6*  (S-1)* 

EETU&N 

EbC 


toaPOiii  P*ucB«a  FcA  FIG  UAL  3.4 

I»Pi.iCil  at.AL*6(A-n,0-L)  , lAiLGLfc  (I-A) 
ccahCA  fc a 1 (1C)  ,FA2  (S) 

LltifcAblG*  FI  (10)  ,o1  (10),U1  (10)  ,V1  (10)  , 

1*  2 (9)  , t>2  (9)  , V2  (9)  ,V(50),T2(9) 
tlliAMAL  *1,P2 
tt*0 
A*  34 

DC  2C  1*1,10 
1=1 

r 1 (1)  =1/Ddw41  (1*3.  141 5 92 65336»7 9323646 2a433) 
CALL  LID  V(Fl,A,2,iAl  (I)  , V , A) 

20  CGMiXMUii 


DC  3u  2*1,9 
1*12  (I) 

12  (l)*DLAP(-i/2) 

CAU.  L1AV  (P2,A,T,FA2  (I),  Y,fl) 
CGMIxaUE 


, . 266 1 9b,  .439004,.  961304,. 
!,.  698462,  .997647,. 992205/ 
,.  *66  329,.  439673,  .561*69, 
!,  .69el3b,. 956133,. 692015/ 


DO  160  AC0PK*1,5 
-Alii  (6  ,4) 

FOAM A2 (Ibl , ////) 

MArtGiA  * 4 t ACOPi 
DO  6 IhAUM  * 5,MA£GIA 
•UK  (c  ,o) 

F0BAA1  (Id  ) 

CD  All  MU  2 

WB12L  (6  ,40) 

FUAHA2  (1U  ,33X,  1 FiGUBE  3.4  CCMFABIbGi*  • , 
1*  Mil 6 Cl  HE  fc  IfiCHAIOULA* ) 

MB  12  A (6  ,48) 

FO*MAl(lU  , 33X  ,45  (*•*)) 

• BI1E  (6,32) 

PO&JU1  tlH0,20X,*F  (X)  ■ 1/6 WBX  (I*FI)  * ) 
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taS II  I,  (0,34) 

b4  FOBtt  41  ( 1h  , 20X,1i  (•-')  ,3X,  • A^PaOXIHAIA  F (1)  *., 

1*  U61BG  BUHABICAL  IBVABSIOM*) 
iUU  (0,  60) 

60  FGB6Al(1tl  ,421,42  {•-•)) 

WB11 A (6  ,64) 

64  F0BHA1(1B  ,44X,17d6IEHFEST*i»  HEX  HOD) 

•BIX  A (6,66) 

06  FOartAX  ( 16  ,42X, 21  ('-')) 

*BATE(6,70)  b 

70  FU«fiAl(lH  , 22  X , 'X  BXAC1’  F (X)  fi**,  12,62, 

1 '6  = 10*  DUEHF.fi*  VEILLoH* ' ) 

»BllE(b  ,72) 

72  FCaaAX(1h  ,21a,'  — *,2(4X, 10  (!-•)), 3(4X,7(*-*)>) 
DO  76  1*1,10 

■fi iXA  (6  ,76)  l,Fl  (I)  ,F  A 1 (1)  ,S1  (I)  ,D1  (I)  , V 1 { I) 

76  F0k6Al(lH  ,21X,lA,bX,A  (HO. 6,41)  ,3  (F7.b,4X) ) 

76  COBI'IHUE 

c 

■BIX  A (6  ,60) 

60  FCiflAI  (1h“/2CX,*  F (X)  * £XP(-£/2)*) 

■fill L (0,62) 

62  kOhib  A I ( 1H  , 2UX  ,1b  ('-'),  10X  ,*  APPfiCXIflATE  F (1)  * , 

1*  USifiG  *UHEa1CAL  IB  VA&olOk  ' ) 

'•fill  A (6, 66) 

ob  rO  fill  AX  ( 1H  , 46X  ,42  (•-*)  ) 

Wall  A (t  ,66)  B 

6b  F Cfi&AX  (Id  ,24X,*I*  ,9X,  *EXACX  F (I)  *,6X, 

I'&XEHFESl  fi**,l2,4X,*  BEaAH Afi*  ValLlOfi*') 

■ BIX  A (0  ,66) 

66  FO«HAX(1rt  , 20a,8(*-*)  ,2  (3X,1b  (•“•))  ,2  (3a,  6 (*-•)) ) 
i/O  94  1*1,9 

■ail  £(6,9  2)  X2  (1)  ,F2  (I ) , F A2  (I)  ,62  (I)  ,V2(i) 

92  F O AH  AX  ( 1 U ,20X  ,F6. 0,2  (3X, Fib.  13)  ,2  (3X,Fd.6)  ) 

94  CobXIBUA 

c » 

■BIX  A (6, 1C  A) 

It  2 FOfifiAX  ( lii- , 2b X,  * *bOX£:  DATA  Foil  btlHCDB  BfifiKAD*, 
1'  BY  Afi  AbTEfilSK  rififiA  XAK  £B  * ) 

•A  ii  A (6, 1C  4) 

1C4  F C fia  AX  ( 1 B , 32X  ,'  FfiCH  A.C.M  AlGOailha  4de,  *, 

1236*  BuaAfilCAL  IfiVEaSlUB  OF) 

•B  IX A (6,106) 

10b  F Oril  AX  ( 1 h , 3iX ,i3 HlAFAACE  XfiAfiSrOfia',  bt  , 

1 *r  AABCGiBA  VAILooB,  CGttaUMI-') 

■B IX  £ (6,106) 

106  fbnrt mX ( In  , 32X,'C  All  OBfi  OF  Xb£  A.C.H.,*, 

1*  Voi. UaE  17,  kUiBAa  10,') 

HfilXA  (6,110) 

11C  FoBHaKIB  , 32X  , • oCiOBEE  1974.*) 
leu  COfi  aIaU  A 
SX'OP 
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nr.  nr 


#h) 

Ilfcfo iei.1  AfcAo*t>4*“H#0-l)  , IMIEGEA  41-A) 

cJhtioA  FAl  4 10)  ,FA2  4$) 

Li hE«6IQa  E(2b),  F (25)  # H(2b)#  V (bO) 

It  (o.tu.*)  GO  Xu  bO 
Ah  * A/2 
Ml)  * 2 
M*)  * 12 
Ml)  * 1 
* (2)  * 1 
40  10  1 * 2,  Ah 

£.(!♦!)  * ( *.  00*  (i*l*1)*B(I)  ) / I 
* 4I*D  * t (I)*I 
10  COaXIJHJ* 
fHh  = Ah 
d(1)  » 2/£  {UH) 

DO  xQ  J - 2,  Ah 

MJ)  * <<J**i'*ii)*MJ))/MAh-J*1) 

20  CGMlIA UE 

ISxGo  « 2*  4 Ad- 2*  (Jlh/2)  ) - 1 
iiO  40  I * 1»  A 

v ( i)  - 0 

of  I ft  Si  * {x  + 1)/2 
JiAbX  * MAO  (I, Ah) 

40  30  J * JFIfiSX,  Jo  Adi 
V(I)  = ¥<!)♦  H (O)/  4f  (I-J*1)  *F  (2*J-IO)  ) 
JO  COAlIA Ui 

¥41)  * IblGM*  ¥ 41) 

16 IGA  * — I6IGA 
40  COMX1AUE 
ft  = A 

50  fA  « 0 

A * • 69 31 4 7 16055* 94b 30 94 1723212145600  / 1 

DO  00  I*  1,  A 

FA  * FA  ♦ ¥ (1) *F 4 I* A) 

60  COMliMUE 
FA  * A*EA 
AETUaA 
EMU 


DGUhiE  FaECISlO*  FUACXIOA  *>1  (S) 
IAFoICIl  AiAi.*6(A-h,0-Z)  , IditUEh (I-A) 
Fl*1/Lb'wA14  6) 

AE'iUAM 
EM  4) 


LOUhoE  FAECI6IOM  FUMcTIOM  ?2<b) 
lflAEICIX  «EAL*84A-h,0-Z)  , IaXEUEA  4I-A) 
F2>4/ (1*2*6) 

AEXUAA 

EMU 


nr  nr 


NMN 


COflBUTLB  t aOO&AU  tOh  t iUU&E  5.5 


IhfLIClX  BEAi.*  8(A-H,K,1,0-Z)  , l a It  via  (I,o,A,ll) 
COafiCA/IAPi>X/EPOC8(10)  ,LAHBLA,  *APSVC,hHu,JToL(b) 
CC«flGA/QUXPUl/vUA£(10)  ,EXPAX  (10,8)  ,IXM10X  (8) 
coaauA/LAPLAC/v  (dO)  , a,  a 
cuaacifc/ahl/osx  Ail  ,S2> 

C0ttaGk/iEHAAP/10*.AbC,lWSUa 

extea aal  pwUAO#paai 

EXXEAAAi.  inAl,Df  5H1,DLFatt1 
h = 0 
A * 34 

C A = 0.0C 

c 

DO  *0  I * 1,10 

CAi.L  LIIIV(PwUAL,A  , EP0C8(I)  ,azix, v , a> 
tUAi)  (I)  = (BHO-1 .)  *APOCU(I)  ♦ ezix 
20  COAXIAUE 

0 

DO  50  J*1,b 

J1CL(J)  = 2*8  ♦ 8 

TGi.ii AC  = 1./ (10.**JTGL  (8)  ) 

CbXABI  = 1 
rxnbua  * C 
C 

DO  40  1=1,10 

CALL  i»l  AV  4BAi.il  ,M  , EB0C8  (I)  ,EZIl,V,B) 

EXP*!  (1,8)  * (EhO-1.)  *EPOCH(I)  ♦ EZiT 

40  CGAi.lAUE 
C 

11  AX  OX  (8)  - IXA5UA 

C 

5C  CoAXImUL 
C 

CALL  IGLoUl 
C 

STOP 

EAE 


BLOCK  LAI A 

iaPLicXS  AEAi.*8  (A-A,K,1,0-Z)  , 1a1EGEB(I,u ,ft ,A) 
CoanCA/lAPUl/£PoCd  (10)  ,lAA8LA,£XP&VC,EiiC,8T0L(8) 
DAI  A cirQCh  /20  • , 4C  . , 60  •,  80  • , 1 00*  , 

1iC0.  ,500.  ,400.  ,800.,  1200./ 

LATA  LAAbu A, El*5  VC,B  80  /I.  00  ,0. 5 5,0.  85/ 

LAO 


SUBBoUllAE  10L0UT 

idi'iiCil  aE AL*b(A-8,K,L,0~Z)  , iAliGEB  (l,«i,H, A) 
COhAOA/lAB  81/EPOC  8 (1 0)  ,i.  AHB*.  A,  LAPS  VC  ,b80  , J1  Ol  (8) 
COttftOfc/O0lPLl/wUAL(10)  ,LAPhX  (10,8)  ,1XA1u1  (8) 
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<xnaoh/A.A*LAc/v  (oG) , a, a 

LLBEaBiod  BESSEL  (3)  , JXOL  (8) 

bEbSELiI)  * 3.9181 
c£SsEL(2)  = b.4ob1 
BESS Ex (3)  = 6.5362 

BESSEL(4)  = 7.4191 
BEbEEi,  (3)  * 6.  1335 

BO  80  *CCPY=1,5 
HAIXE  (6,4) 

4 FOkBAX  (Idl,////) 
bAbiil*  * 4 ♦ a COPY 
i)0  8 xBBGii  - S|ttA2UI< 

.£  11  £ (b,6) 

O iQhbkl  (Iti  ) 

8 COM!  HUE 

• 6 11'  £ (6, 1*) 

12  FOBBAX(1a  ,23X,’EIOU&E  3.5  EFFE CX  of  8a«ioa-kAP8S CA 
1*  SEAFtth  •iOiEBAbCE*) 

ifill'c  (6,1 4) 

14  FGEflA5L(lB  ,29X, 53  (•-*)) 

■ a I XL  (6,  1 o) 

16  FOafiAT  (180,  1C*X,37HltiIS  TABLE  CuaPAdES  CGaEBAX'S  BE  SSaL  , 
1 1 F OaClxOa  aEbULTB  Al  FIVE  EPOCHS  mITH  LaPaACE*) 

m a lx  L (6,18) 

1b  FGBbAX  (1H  , 10X,'1MV£aSI08  EBSU..XB  FOB  Bulb  X8a  EXACl  • , 
I'WliaLkAXAC  SOLUX1CA  OF  Xd£  FUBCX10MAL  A*L  1 BE  ') 

■ LITE  (6,20) 

Xb  FOE MAX  (Id  , 10X , 4 APFaOXlM  AT E M EaXO*-A AP HSofc  SEAdC d », 
1'SOLUXi.GB  OF  THE  F U NOTIONAL  Nlid  VAeIOUj  XOaEiANCES.  •) 
W&XTE  (6  ,30) 

30  FOaBAl (InO, 10X,*d/a/1  BEAM  SEaVEE  LOAD  , ' , 

1*  LABBDA  * 1.0,  E (B)  - 0.9b,  4.  - 0 . • ) 

DO  BC  JJSIkSl*1,S,4 
sLAST=OFIkSX*3 

■BIX E (6, 3b)  a 

38  F Oak AX  (1dO,3oX,'XN  VEkSlON  OF  LAPLACE  X kafcSF  CFua,  • , 

1228  SXEdFESI'S  HE 1800,  N*,12) 

■kITE  (6  ,40) 

40  FOaBAX  (la  ,1BX,'S0BB  OF  *#73('-')) 

■k  IX £ (6,42) 

42  F 0BhAX(l8  , 19X,*  BESSti.  yUauBAXLC  XOaEaANCE  FOa  », 

1 •kENXCN-aAPHSON  SEAaCH  SOLUixON  OF  FU  NCxiOJi  AL' ) 
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• H1XE  (0,44) 

44  FOaHAI  ( 1h  , 10X,  'EFOCH  FUHCXI06S  5JLUXI06  OF 


EO  4 b o =JF i&SX  , Ji.ASX 
46  C06T16UE 


•KITE  (6,50)  (JlOi.  (J)  ,J*JFIaS1‘,JLA3X) 

50  FOuilAl(lH  , 1C  a , ' X (COLEbAH)  FU6CX106  aL 
14(*10 **-•  ,11,61)  ) 


4 A lit  (o  ,56) 

56  FuhdAl(1h  ,ICI 


• aIX E (c,  oC)  EPOCd  (I)  ,3 ESS  EL  (I)  ,*UAl>  (I)  , 

1 (EXPmX  (I,J)  ,j=JF1 HSX,JaAST) 

FCKMAX(1d  , 1CX  ,F5.0 , 5a  ,F 7.  4,  5 ( IX  ,F14 .1 1)  ) 
CoKllKUE 


it  HUE  (6,  64)  EFGCri  (I)  ,Q  UAD  (I)  , (EXFiX  (I,  J)  , J =JF  IASI,  JL  LSI) 
FOBhAT  (1H  , 10X,F5.C,10X,5  (1a,F  14.11)  ) 


•BIX  E (6,70) 

70  FGBd  AX  (1h  , 10a  , 'TOTAL  6U06EM  OF  AEUXON-HAF US06 ') 


• All  x.  (6,72)  (I  Id  Xu  X («i ) , J= JFIaeX  , JLASX  ) 

7 2 FGBEAX  (Id  , 1GX  , ' 1 1EB ATIUMS  TO  EVALUATE  10  EF-OCdS 
14  (14,1  IX)) 


6C  Co  61 1 HUE 
aEXUBK 
EKE 


MJaaOuxIeE  LIKV  (F  , 6,  T, FA,  V , a) 

IM ELICIT  EFAa*8  (A-b.u-Z)  , IKlEGE* (I-K) 
C0rtbO6/LArL  AC/EFG  Ch ( 10)  ,AFFaoX  (1C,o) 

El  ME  waI CK  E (25)  , A (25)  , 3(25),  ¥(50)* 

IF  (xi.Ew*)  00  To  50 

Kb  = 6/2 

E(1)  = 2 

fc<2)  » 12 

MD  * 1 
M2)  =1 
uu  1 0 I - 2,  Mb 

a(I*1)  * ( 2. 0 0*  (2*I»1 ) *E  (I)  ) / I 
Ml*1)  * F ( 1)  • I 


1 


gags 

■" - - 

■ 1 

10 

Cu4  J.  14  UK 

P4fa  * fed 
u(1>  = 2/e  (lid) 

DC  20  J s 4>  | ti  h 

b(J)  * IMh-J ♦ 1) 

20 

COH'iidOfc 

ISXGft  = 2*i*d-2*  (4*i/2)  >-1 

00  40  I * 1,  4 

V(I)  * 0 

JiluSl  = (i*1)/2 
obAd*  = dlk  0 (I  «k  d) 
tO  30  J : oHfibi',  Ol  ASi 

V(I)  = V(I)  ♦ d(J)/if  (I-J*1)  *F  t2*J-lO)i 

j 

30 

CuftXliUE 

V(I)  = idled*  V (I; 

XdlGft  = -iSIUM 

40 

comiaui 
n = * 

50 

Jr  A * C 

A = . o»3l471d0bb9&4d30S4l723x:1214:>8£>C  / 1 

to 

LU  60  i = 1,  4 
k A = k A * V <i)*P  (I*A) 

LOH'x’ 

' I 

r A * A*tA 
adXUdft 

;| 

c 

d4b 

c 

kA,Uo*.k  trade  IbXCM  Fbfteiluft  Pwdab(d) 

ittPi.itIA  *£Ae*6  (H-fc,K,  1,0-2)  , Iftloeoa  ( l,d  , A , ft) 

Codfieft/iftPbl/EPOCd  (10)  ,LAfld*A,EiPdVC#dtte,Ji.ot.(tf) 
o=1-  C.XPSVC*  (e  AHbbA  + S) 

eMdG  A = (-b*bSi)4T  (b*b«-4*S*EiP.»t#C)  ) /(2*EAt'bVC) 

luUAb*bEXP  (-GB  EUA*a)  / (S*OttiU»A) 

aElUdft 

Eli  0 

c 

DOUdeb  PdiCidIG*  HJ4CTI04  Pddl(b) 

IHPtiC  XX  fat  A a.  * fa  (A*ii|K/Lf  ; Xft  X oGEd  (if  Of  A/  4) 

COftdoA/INt  u'X/EPoCd  (10)  #LAdbbA,  dit-d  VC  ,ihu  # JAG 2 (o) 
eGflfiCh/ftal  /OdA  AJAX  #SS 

EXAEMhaL  iCKt»1,DFan1 
dd  = d 

CAifi.  Zi.de  (FOdft  1 , b Arid  1 , GSXAaT,  GHdG  A) 

OdTAEX  = OtiEGA 

?«Bl=DfcA*(-C«EGA*Z)/  (S*JttdGA) 
bEXUad 

END 

c 

c 

LOUded  PadC  loICft  PUUCTiOM  Fddl  (d) 

iHtoieiX  adli,*6  (A-H,  K,  1,0-Z)  , i«ldGE4(I,J,a#k) 

CGBdGfc/IftPbl/dPGGd  (10)  ,oAflBiM  , EX  Pd  VC  , dd  C,  o *Oo  ( 8) 
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on  r;  r r r>  r. 


f hoi  * b - XAHBDA*  (1- (V(1*i>*tXP;»iJC)  )) 

EEXU  Hk 

EHD 


DC-Ud^t  f HECIdlQH  fUHCl'lOH  DP flat  (bj 

IdPliCiT  HEA  i.*6(A-H,K,E,0-£)  , IHIEGEfi  (1,0,  II,  k) 

COhdOH/lkP UT/EPoC H (10)  ,LAHB0A,  EXPSVC,ktio,JiOl(8) 

bftthl  = 1 - i.AdBDA*EXPbVC*  (V  (1*S*EXPSVC))  **2 

HEiO.uk 

Ik  0 


DGUdiE  ifiEClblCH  f UkCliOk  dOiHhl  (d) 

IhPi-lCIT  *EAL*b(A-H,*,  L,  0-2)  , lalEGEE  (I  ,G,h,k) 
Co«h0ti/I»PUVEPuch(10)  ,LKHbDk,  dA  kb  VC , HbU,  OXOE  (b) 
vUt  Ha  1 = 2*LAfl£UA*£lCPSVC**H*(1/(Ud*i.APdVCj>**3 
aEI'UHH 
EH  £ 

DO  U bib  r Hi  CIS  I OH  PUHCTIuH  fUah  1 (OEEGA) 

1C111  aEAL*  b (A  — d,  K,  L,  0~X)  $ IklEG  EH  (X  ,0, a , H) 
COhdCk/dal/GdT AkX  ,bS 
iOMEl  = t hilUOrtEO  A)  - Sd 
HLi'GnH 
EH  L 


dUddcUilHE  ZcblO  (t  ,DP  ,X  bXAHl,  aXEhG) 

xEEaxcEI  HE  Ai«* b (A* h/ Kf  E| 0*E j f IHiEifEH(I,d,a,k) 

COBhCH/kEHk  Ac/IuaHhC  ,IlkSUfi 

IIEHaT  = 0 
X = XSl'Ahl' 

1C  COHilHOE 

XktH  - a - f(X)/«i(X) 
aELAcC  = (aNEH-X)  /XNEu 
a = aHE ■ 

I'XEh  AX  = ii  EH  AT  ♦ 1 
If  (IXEiAX  «G£.  bO)  GO  XO  5U 
if'  (LAdd  (HELaCC)  ,a.E.  iOLHHc;  GO  XO  20 
Go  io  10 
2C  Coax la Ui 

iiHSOh  » i'iuaoa  ♦ iieeax 

XEEhc  = X 
aEXUEH 

50  COHaIHUA 

«£  IX  E (b,bbj  HEIacC,  lOXHHC 

d5  fGkdAX  (Id  , *ITEHAX*bO,  kELACC*'  , 1)24 • 16  , ' , ioLHkC*  • ,020.  1 b) 
IlkdUfi  a IXHSOrt  ♦ IXEaAl 
XXEfiO  a X 
HEXUaH 
EH  L 
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COflEUliM  FaOoMAE  Eo«  iiUUME  J.  A 


1 


IMELICII  tiEAL*8  (A-fc,o-Z)  , IMiEUEM(I-li) 
cobttcti  /Aiioz/  MuO,*. 

LlflEASlO*  SEaIXI  (11, 8)  ,SE«IX2(11,d)  ,V  (50),*SIM  (8)  , Muuli,  (11) 
EiXEMMAi.  P^UMLl,  FwUAL2 
L Ai A ZS2  A / «2|t4|«0|t6,1*|2i,J«,4«/ 

LATA  bu olAi  /«b|  «o  , • 7 , • b , • 9 ,1*  1 ,1*  2, 1*3, 1*4,  1 *0^ 2*/ 

t 

h=0 

M=34 

C 

00  U IfihOs1, 1 1 
MHO  = MdGl  A (IMHO) 

AXFCXm  = LAMS  (1.- MHO) /MHO 

Do  12  JZ=1,8 

Z « ZSI  A (JZ  ) / Alt  Clfc 

CAx.E  JL1MV  (FyLALl , A,Z,£Z1Z,  V,fl) 

SEt.111  Undo  ,UZ)  - E*1Z*MTFCTM 
CkLu  LIMV  (j?WUAC2,M,Z,KZIZ ,V,fl) 

SEZI1  2 (I*nO,JZ)  = EZIZAMTFCX* 

U COMTIMUE 

c 

LO  Mb  MCoE¥=1,5 
MM  IX  L{t  ,20) 

20  r OMfl A1 ( lb 1 ,////) 

BAMUlM  = 4 ♦ t»CO  f i 
LO  aO  I MHO  A = 5,  fl  AfiOld 
■fi  IX E (A  , 24) 

24  FGfih  A'l  (Iti  ) 

28  COMXIMUE 
C 

AM  IX  E (6,32) 

32  t o Md  AX  ( 1 d , 2bX F I CUBE  3. 6 CGMmECIIGMS* , 

1'  iCn  LISCoA'ilA  UO  US  ilMSX  LEMIVAIIVE1) 

■MIXE  (6,36) 

3b  FGMriAl  ( 1H  ,2 ol, So  (*-•)) 

• BITE  (A ,40)  A 

40  FGMiUX (1b-, 33*,* AFFMGXIfl Allots  USInG*, 

122b  SXEHFZ.al'S  ttElaO L,  A*, 12) 

AMIxE  (A  ,44) 

44  FGBbAX  ( 1H-,  20a, ' &/A/1  SC  A a.  EL  BEAM  CUflUA.AH  VE  • , 

141  h IDoiMESS  AX  X 1 *AMS  (HU)  *Z*  (l.E.  , AT  1-Z)) 

WM1XE  (A  ,4b) 

4b  F0Bb*X (IdO, 47 X ,*  mI'XHOUX  COMMECIICm  XEBfl*) 

■ SHE  (A, 52) 

52  AO  EH  AT  ( 1 dO  , 45  a , * SCALEL  IMITIAi.  S EM  V EM  LoAL*) 

AMIXE  (A, 56) 

56  fOttflAX  (Id  , 19X,7b  (•-•)  ) 

wmIXE(A,o6)  (ZSlM(JZ),  JZs1,b) 

60  FOmHAMIH  , 15X  , 'a  bo  * ,5  a,A  (F3.  1 , 71)  ) 

■ MIX  E (A, 64) 

64  FGMAAX  (1H  , 15*  b (2X,  8 ('-«)  )) 
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— 





do  it  iaao*l,ii 

W&IXMCfBO)  fiHOIH  (1680)  , (SbZI'ii  U«80,Jo)  , J**1,8) 
Ob  iObbA2(lH  ,l5X,Jr8.  1,6(2i,F8.o)  ) 

12  cOMXiBUJS  j 

c 

kkWL(b,lb) 
lb  £0B6Al(1b  ) 

•BIT i.  (0,44) 

66 IT t (0,80) 

80  Jr  OBb  A2  (IhG  ,482  ,'  «i  1T  b C0HMECII06  llSEfi*) 

W6I2E  (6,62) 

UfellJL  (6,56) 

bitl££(0,oO)  (ZSIB(JZ),  JZ*1 , 8) 

MBIT £(0,64) 

BO  84  1680-1, 1 1 

*6  IT  £ (6 ,68)  6 HOI  6 (1680)  , (S£oll2  (16  80,0*)  , J2«1,a) 
64  c062l6Ui 
6b  CObl'lMUE 
Si  OP 
16  D 
C 
C 

SBbBOUX  IBS  lDiV(P,o,l,PA,V,Q) 

IMPLICIT  UoAjl*  6 (A-H,0“Z)  , l6TJsOEt  (1-6) 
COfid06/*APi.AC/EPoUi(10)  ,APBhOA.  (1C,  6) 

Bidi.*i>10«  *(25),  6(25),  H(25)  , V (50) 
it  (d.bwft)  00  10  50 
Bh  = 6/2 

Ml)  * 2 

£(4)  * 1l 

Ml)  * 1 
M2)  =1 
00  10  1 - 2,  Mil 

& (1+ 1)  = ( 2«i)0*  (2*1*1)  *£(1)  ) / I 
1 1 1*1)  - k (I)*I 
1C  006 216 UE 
p«u  - 6b 
8(1)  * 2/1(68) 

B0  40J  = 2,68 

8(0)  * ((J**P6i)*£(J))/*  («8-J*1) 

40  00  62  168* 

IS  106  * 2*  (6b- 2*  (68/2)  )-1 
00  40  1 * 1 , 4 
v(i)  * c 

Jl  1632  * (1*1) /4 
02 AS2  * dl60(I,68) 

BO  80  8 * JflAS'i,  J2AS1 

MI)  * V(l)  ♦ H(J)/(f  (I-  J*  1)  *1  (2*J— 1*1) ) 

80  COM16U6 

V(l)  * ISIUM*V  (1) 

1*106  - *18106 
40  006X1608 

8*6 
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»A  * 0 

A = . oS  3l471b0io^i»M5i0V4  / % 

IX)  i»0  I*  1,  n 

*A  * iA  ♦ V (iJ*E  (1*A) 

CQtili.AU  A 
Jr  A = A**  A 
hElUak 
ch  C 


UJj OLE  daACIblOtl  JUNCIIOA  *wUAUl  (S) 
IflEi.lCI5  «EAa.*6  (A-b#o-Z)  , 2.N1EQ&E (I-») 
CCiHflC*  /&H02/  Edo  ,4 
b = 1 - BdG  - S/A • 

ONEGA  = -d*lL>uBT  {d*d*2*)») 

tgUAil  - LEXt  (-OdEGA*A)/  (S*o£iEUA) 

BE1USE 

END 


DOUUi.E  PEECISIO*  JUJICIIOA  BWUAU2  (S) 

IhBLiCIl  *EAL*0(*-H,Q-2)  , IaXEQEa(I-A) 

COttauA  /EriOZ/  EdG,A 
b = 1 - Edo  - S/2. 

OhtGA  * -b»DSU*l(B*B4-2»S) 

f'^UAi)2=I»£,X£  (-GHEG  A*Z )/  (d*3d  &GA)  - LfcXP  (-**  (2.*kdO*S)  )/  1S*S) 

aETUaA 

EMC 


! 
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r.  c 


cOMPUTEE  PaOG*AM  FOE  APPENDIA  TAELlS,  ktiO  EOT  x,uUAL  TO  OHE 


IMPLICIT  ELAL*8(A-h,K,L,0-L)  , I* SEGEE  (1,0,  M,M) 

COMMON  /INuUl/  1 V A*.  (1 6)  ,K  V AL (7)  , 6C*» HI  (10  ,1b  , 3)  ,EhOlE, 
1LSIE  (3) 

COMMON  /L  At  LAC/  V (50)  ,h,N 
COMMON  /NEHEAP/  i'QLEEC,MAi.ITN 

M*C 
E*  34 

MAXITN=20 
lOLEEt*1.  D-20 

DO  bO  I PAG  6*1  , 3 

HEAD(b,2C)  EHOIN,  (ESIN  (IIASlE)  , ITADLE*1,3) 

EEAD(b,2C)  (TVAl(O)  , J*1,b) 

*EAD(b,20)  (IVAL(O),  0*7,12) 
aEAD  (b,  20)  (IVAL(O),  J=1 3, 1 6) 

FO Ah  AT (6  * 10.0) 

DC  50  ITAELE=1,3 
CALL  ELWNE  (IT  ABLE) 

CALL  ElMLI  (ITABLE) 

CALL  aEMEK  (ITABLE) 

CALL  EL  GAM  (ITABlE) 

CONTINUE 
CALL  n20Ul 
CONTINUE 
STOP 
AN  D 


bi.CC K DATA 

IMPLICIT  BEAL*  8 (a-h,  K,  L,  O—Z)  , INTEGEE  (I,0,tt,N) 
cOMMofc  /INOUi/  T V AL ( 16)  , KV  AL  (7  ) , SCENT  (10,16,3)  ,hdOI N, 
IlSIN  (3) 

LATA  KVAL/4. ,3. ,2*, 1.5,1 .,.5,. 2/ 

EEL 


SUBEUUTIeE  LINV  (F,h,X,FA, V,M) 

IMrLlClT  «EA^*6(A-E,0-L)  , INTEGEB(I-N) 
COMMCN/LAPL  AC/ EPOCH  (10)  ,APPaG*  (10,6) 
olEEASlGN  E (25) , F (2b)  , H(2b),  V(50) 

IF  (O.Ew.E)  GO  TO  50 
Ah  * e/2 

MD  * d 
E (2)  * 12 
Ml)  = 1 
F (2)  * 1 
lC  10  I * 4,  EH 

E ( 1*1)  * ( 2» uO*  (2*1*1 ) *£  (I)  ) / I 
F(i>1)  a F ( 1)  * I 
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r r c.  r> 


1 


10  COMliMUE 
PUb  * lid 
H(1)  * 2/F(ka) 

DO  aO  o 3 2,  kH 

d(J)  * UJ**PH&)*£(J))/F(MU-o*1) 

20  CGkxIMUE 

ISXUM  3 2*(kd-2*(kd/2)  )-1 
u(i  4C  I - 1,  I 
V(I)  * 0 

Jilka't  3 (l*1)/2 

JLASa  * AlkG(I#Md) 

DO  30  o 3 JFlaSl,  ooAil 
V(I)  3 V(I)  ♦ d(J)/(F(I-J*1)  *i,(2*J-I*1)) 
30  COk'IIkUE 

V<I)  * ISIGli*  V (1) 
iillik  3 -ISIGM 
40  OlkXlkliE 
a * k 

AC  FA  * 0 

A * . 63J1471«055i»94A30!#4l723^l2l4AttD0  / 1 

00  60  Is  1 , k 

FA  * FA  ♦ V<I)*P(I*A) 

60  CGMXIkUA 
FA  * A*FA 
££1U£M 
Ek  0 


EUbfioOilkE  ttakkfi  (IXAB££) 

IMPoICIX  it£Ai.*6  (A-fc,  K,i#0-Z)  , Ik'j.EGEk(I,,j#h,b) 

COHHOM  /I MO  Li/  X VAi.  (16)  ,K  VA- (7 ) ,6c  Lk*X  ( 10,  1 o,  3)  , IihGIk , 
1261k  (3) 

CO  AH  Ok  /(.ArlAC/  k<60),M,N 
COHflCk  /AdOZO/  AdG,Z,06TAET 
EAZEkkAi.  Pdka 

6CALIbG  fell M dlGagb  * 1 
kUO  3 knolk 

Ai.Fi  gd  * (1-kHG)  • (1-MHO) 
klFClk  3 DA05  <1-*dC) 

2 3 kilk  (II*61£) /UZFCXd 

00  *C  o*1#16 

1 * iVAi.  (J) /A LF 6 wM 
CAL*.  i.2MV  (PWk£,k,7,T££M3,V,H) 

£2 MX  « 2 ♦ (£U0~1  •)  *T  ♦ IEEd3 
SC1HX  (1,J,1TA6L£)  « £2  MX*  MX  fold 

aO  Co  k T lk  U E 
C 

feEZUMM 

AMD 
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c o r r.  r r r n r 


DOUBLE  FiEcIilGM  FUMCXIDM 

IKpLICll  EfiAi.*6(A-H,K,I.,0-z)  , IEXiGEB(I/J,fl,E) 
COBBCE  /EhO 20/  EUC,Z,0SXAEJ 

ChiU A * -(1-EHO)  ♦EiyEX  ((1-ihu)*4l-ihC)  *E*a) 

A'W HM  s DEEP  (-OBEG  A*Z)/ (S*OHLUA) 

hi 1U*E 

EEC 


DUBEuUliMi  BZhDl  (1XABLE) 

IhPLICIX  a*Ai*8  (A-H/K,L/0-Z)  , IhlLGEa(I, J , H ,E) 

CGhflOE  /I ho  bl/  X VAi.  (16)  #KVAi.  (7)  /SCj.M X (10/  1 0,  3)  , BHOIE  , 
IZSlE(J) 

CUhhOh  /LAPLAC/  V(SO),B/M 
CGhhCE  /EHO  ZQ/  BUG  / Z /GBXABX 
iXXEaEAi.  PBDl  / FMD 1 / DFB  Cl  # DDFhD  1 

BC  Aj«1EG  MIX'  fa  BlGS^E  = BUG 
EEC  = BHOlh 

mLFSwE  * (1-EHO)  * (1-fittQ)  /EHQ 
EXPCXi  * DA6B  (1-HHO)/EHO 
OSXAEI  * 1 

i = Zalh  (ll'AFLE)  /WTFCXE 

LuZ  = DiXP  (-*HO*Z) 

IF  (EUG  .LE.  1.C)  GO  10  10 

CAIx  ZLBO  (0FBD1/  GDFB  D1  / 1 D-  10  /SBAk) 

OELXA  = -FBDl(SbAE) 

ilAHl  = Hb AB  ♦ DaIXA 

CAii  ZLEG(F  BOl/DFBDl/STABX/OBXABX) 

10  GOhXlEUE 

DO  10  0=1/16 

X * XVAi.(J)  /ALFSjB 

IF  (X  i L£«  ZJ  EZx  X = 0 • 

IF  (1  .01.  Z)  CAoi  LIE V l tflol  / E/X'/UT/V/h) 

IF  (1  • 0l«  Z)  ALII  = HI  ♦ (X-Z)*Ei.Z 
LZ«X  * Z ♦ (EU0-1.)*X'  ♦ EZIX 
BCLEl  (2/J/IXABLE)  = EZHT**TFCT* 

10  COhXIEUE 
eE'XU  EM 
EMD 
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nr  n n nr  nn 


. 

U)Ui>*.t  PaLiXSIGA  t UMCTiOM  PaEl  (S) 

AEJL.*8  (A-0,K,  i,0-ZJ  , jLAlLGE*(I,d,A,Aj 
CO AAl*  /aduZo/  Ado,Z,OalAAl 
CCAAcM  SS 

LX'XaaA*!  FOAL  1 fl)t  AEl 
SS*  S 

CALi.  Z**0  (*  OAOl,t»PA  Li, GSTASX, OMEGA) 
oSlAkl  * CAECA 

PALI  * U^At1  (-GALi»A*Z)/ (S*OALUA)  - EAXP  (-Z*  (*AO«S)  )/(S*») 

AE1UAA 

EkD 


COUBA.E  *AECiaIGA  FUMCIIOM  FALl  (a) 

IMPIaCI'A*  *EAZ*b(A-d,K,E,0-Z)  , XA1  AGE  A (1,0,  A,  A) 

CO A AC A /fiHCZu/  AaO,Z,OSiA£T 

FA  1)1  * S-kttO*  (1-CEXP  (-S)  ) 

aEXUAA 

LAD 


DOUBLE  PaEClSIGA  EUACXIGA  i)FADl(S) 

iA^LICII  fcLAL*8(A-b,K,l,0-Z)  , laitGEB  (I , J,  11,  A) 

COAACA  /aAOZO/  kd  G,  Z , 0 ST Akl' 

LFMDl  * 1-«bC*EEiF  (-S) 

BEX USA 
LAE 


EOUbaE  PkACISIOA  FUACXIOM  DLFADl(S) 

1APZICI1  kLAZ*b(A-ta,K,i.,0-Z)  , lAlEGEB  (i,u.  A,  A) 

COAACA  /aaozo/  mhg,z,ostamt 

CDifiCl  « AdO*L£XP  (-S) 

kElUAA 

LAE 


COUdZE  PkLclSIOA  FUACTIOM  FoaLl  (OMEGA) 

IHrE^CIT  AEAi.*8(A-8,K,L,0-Z)  , IAILGLB  (1,0,  b,  *) 
COAACA  aS 

FOAE1  * <Az1(OAAGA)  - AS 

AS  1 USA 

AMD 
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r'  r.  C r.  n r n 


C 


aGbhCUXIliZ  jiZhEK  (I1ABE  £) 

IhPEICIl  EEAL*8(A-h,K,E,0-Z)  , I*lEGE*(i ,u,h, ji) 

C CtthoE  /I«GLX/  X V Ai.  ( lo)  » KV  AE  (7 ) , ECEal  (10,16,3)  , A HO  la, 
iZElb  (3) 

CGtihCA  /-Ar LAC/  V(bO),fl,N 

cooacK  /ahozok/  hbG,z,GEXAai;#* 

EXXE  aJUL  PHEK,  FH  *K  ,CFH  EK  ,DDPhEK 

EcAEIEG  U IX £ ElGSufc  * (K*1)*ahO/K 
EhC  * EhOI* 

EG  30  1*3,  b 
K * KVAL  (1-2) 

AEFEwh  * K*  (1-KbOJ*  (1-aaO)/(lA*1)  *iUO) 

HIFCSE  * K*DAbS(1-EhO>/(  (K+1)*fthC) 

uSlArfl  * 1 

Z = Z3Ifc  (IiABxx.)/«XECXil 

EEZ  = DEAP  (-Ah 0**.) 

if  (AUG  .EE.  1.0)  GO  SO  10 
CAi.E  ZE*0  (EPflEK,GEiHEK,0,SbAA) 

GfiLXA  * -EhEK(ShAfi) 

El  Ail  = ShAa  ♦ OEElA 

CAEE  ZEEG  (PfltK,!,*  HEA  , El  A E1‘,  OEiAAI) 

10  GOKXIAUE 

DO  20  J=1,16 

X * X VAE  (J)  /A  LFSwh 

IF  (X  .EE.  Z)  EZ..I  * 0. 

IF  (X  .GT.  Z)  CAL-  LIE  V (PhEK,*,X  ,hT  , V,  h) 

IF  (X  .GX.  Z)  EZil  * hi  ♦ (T-Z)*E EZ 

EZbl  * Z ♦ (aHC- 1 • )*'£  ♦ EZII 

ECEHX  (I,  J , X ZABLE)  = £ZWX*UXFCZfi 
xO  CONTINUE 
30  GukllAOE 
EE 1 UAH 
EhO 


LGUe-z  PaEcIEXO*  FUaCAlOA  PhEK  (a) 

XflPzIClX  &EAE*6(A-h,K,L,0-Z)  , lalEUEfi  (I,  J,  ft,  M) 

COHhoA  /atiCZGK/  ahO,  Z,UEl’ARX  ,K 
COKftOh  EE 

EllEafcaE  FuflEK ,DF  fl£K 
SS*  E 

CA  Li.  Zz*G  (F  ChEK,  UFhE K,05l  Ail  ,Oh£GA) 

GElAE'X  = GhEGA 

PAEK  * LEAP  <-OfiEGA*Z)/ (E*OHEGA)  - DEEP  (-Z*  (ahO *E)  ) / (E*e) 

«E1UaA 

END 
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00UO4.C.  *£aCj.E1Ga  * JaCTIGN  (BU 

l*k+±Cll  a£AL*64A-H,K,  L,0-I)  , IA1EGE*  (I  ,o#a,A) 

cOMttuN  /AbuICK/  ufcC, £,CSlAtt,  A 

FHaA  * S-AuO*  )**K) 

aETUaN 

LA  L 


UGUaLE  AMEC  Li lOti  FUNCTION  Ui't U.M(S) 

LkiLU.il  kfci*.*8  4A-fa,A#i,C-£)  , IATEGEN4I,J,A ,N) 
CGNftcA  /MhGECK/  A UG  , E,CSXAJll  , K 
ut££*  « 1-aMO*  (A/  (K+S)  ) **  (AO) 

«£XUMA 

EaL 


DO UbLL  PRECISION  FUNCTION  DuENE*  (£) 

lkkt.xi.li  ££*1*8  (A-t,A#i.#G-Zj  , XAXaG£E(I#o,£i,A) 

ConttcN  /AhOZCK/  £ rid  , Z,<JkTk*l  , K 

DOME*  * hhC*  (KO)*  (A/ (£♦£)>•*  (K>2)/K 

aETUEM 

Ui  L 


DOUtict  ItLClolOA  FUNCTION  fcOAEA  (Cb£UA) 

IktxLCLl  N£JU.*b4A-h,K,  L,G-Z)  , iNlEGfih  (I  ,ii  ,N) 
CCttflcA  £8 

JroBEA  * £ b£K  (GH£G  A)  - SS 

aETUaN 

END 


aUBaoUXIAE  BEGAN  (IXABLE) 

mEAL*8  (a-B,K,I,0-Z)  , iMXEGEM(I,J,b,A) 

COMMON  /INGUX/  T VAL  (16)  #K  VAL  (7)  ,£C  L«S  ( 10,  1 8,  3)  , Ad  Ola  , 
1 ASIA  «3) 

CuNfiGA  /LA*  I AC/  V (5 0),tt,N 
CGHnCN  /MfluZO/  AdO,Z,OSTABT 
EXTaANAI  iG  Aft,  J?G  AH , DFG  AN 

SCALING  All  fa  AIGSwN  * BUG 
AriC  « EtiOIA 

ALFNwa  * (1  -iUd)  * (1-EHO)  /Nto 
•TfCXM  * UAoE  4l-AilO)/AU0 
GaTAaT  * 1 

A » *£IN  (IT  ABLE)  /HTFCTa 


c 


Iff  (aUC  .La.  1.0)  GO  10  10 
EbAk  = Etto  - 1 
DELIA  * -FGAh  (SLAB) 

El  A a'1  * SLAB  ♦ DELIA 
CALL  2E10(EGAh,DirGAI!,BlAkI,0£rA£l) 

10  COEklhUE 

DO  20  J*1 # 16 
1 = 1VAL(J)  /LLiSutt 
if  (i  . LE.  i)  Ei.il  « 0. 

IF  (1  . Gi.  2)  CALL  LIBV(*GAa,B,X,EZIl,V,h) 

E2«l  * 2 ♦ (BUO-1.)*!  ♦ E 211 
SCLBk  (10#J,riABi.E)  = EZHl*  Ki'FC'fB 
20  COBlIkUE 
C 

nETUBE 
E AD 
C 

c 

DOUtu.c,  tkZCISIOB  JrUMC'IIOM  PGA  h (E ) 

lapiiwia;  AEAL*8(A-b,K/L#o-z)  , iaiegea (i,j,a, e) 

COHtiOM  /Bho20/  AtiC/if OS'IAfiT 
COtiflob  EE 

LA  Xh  EttAL  FOGAM  #i»r  GAM 

£>E  = E 

CALL  LEAD  (E  UGAB,D  rGAfl,  OSIAkl, OMEGA) 

CEIAa'I  * Oh  IGA 

PG Ac.  * LEAP  (-OMEGA* Z)/ (E*OHhGA) 

aEIUaE 

EE  D 


DOUELE  PBECISICk  fb* CTIOd  FGAM(E) 

IMPLICIT  BEAL*8(A-U#K,L,0-Z)  , aEIEGEB  (I  ,J#tt,k) 
COtittOE  /BbOZO/  tobC,Z,OSTABl 
EGAh  * S-ahO*DLuG  (1*S) 
d El'll  EE 
cJiL 
C 
c 

DOULaE  PBECIeICE  EUECTIOE  DFGAB(E) 
ittPLICIX  LEAL*  8 (A-h,  K#  L#  O-Z)  # Idle, GEE  (I,u,h,k) 
COflhOE  /AbO  20/  AbG#2|0aI'AJi7 
UFGAa  * 1 -BHG/  (1  ♦ S) 

HEIUaE 
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UOllo*.E  EaEC lEIOM  FUNCTION  FCGAd (GdEG A) 

IMaeaCIX  XEAL*6(  A-d,K,L,0“Z)  , IaXEGeN  (I,u,  a,  *) 
cOBBCN  /AtiOZU/  «nC,E,ObXAfil 
CuBduN  56 

FOGAM  = FGAd  (OMEGA)  - S6 

«ETU  NN 

END 


EUBEcUIIeE  ZENO  ( * , DF  ,X5XAeX  ,XeENO) 

IMr^xCI'i  EEAE*8(A-h,K,L,G-Z)  , INIeGEB  (1,-U,  6,  N) 
COMMON  /NEKAA*/  T Gi.fi NC , d AX1 XN 
lTEaAl  * 0 

X = 151 A dT 
10  CONTINUE 

AM E ■ * X - A (X)/DF  (X) 

BEeACC  * (XNEM-X)  /XNE* 

X = XNEa 

IX EE AX  = 1a  EE AT  ♦ 1 
IE  (I2EBA1  .GE.  dAXIXN)  GO  XU  50 
IF  (uAub  (BEEACC)  .IE.  xOLBNC)  GO  XO  20 
GO  1C  10 
20  CONI  IX  UE 
XI EEC  = X 
EEIUeN 


50  mE IX E (6,60) 

6C  FOEEA1  (Idl  , * N E»IO*- BAFhSGN  oEAECd  IN  EUoNwUTINE  ZENO*  , 
1*  EXCEEDED  MAXIMUM  ALLOWABLE  NUdBEE  OF  aXEEAXIONj. •) 

WE IX E (to, 62)  MAXI IN, EEL ACC, XOEaNC 
62  FOFBAa  (IbO  ,*MAXIXN  * * , X4,///,  • eEEACC  = »,0<.8.1o, 

1///# ’XULeNC  = *,D2b.l6) 

SlwP 

END 


bUto&UUllNE  oeuUX 

INEe ICIT  aEAE*6(A-h,K,E,  O-Z)  # IXTEGEE  (I  ,u, 6#  N) 

COMMON  /IN  OUT/  X V AI  (16)  #KV  Ae  (7 ) # 5CEM1  (10,16,3)  ,NhOlN, 
IZblN  (3) 

DIMENSION  ISCLNX  (10,16,3) 

DO  »C  NCOfX*1,5 

NK HE  <6,10) 

10  FOBM  AX ( 1d1 ,/// ) 

BA  EG  IN  * 4 ♦ NCOP  1 
DO  1a  IdltG  M *5,  MAN  GIN 
NB I1E  (6,11) 

11  FONMil(1a  ) 

12  CONTINUE 
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UO  60  ilAbLt=  1/3 


SCALCB  * 1000. 

Tt  (6Ct.N1  (1,l6,IX  ABLE)  »G£  • 10.)  SCALOE  * 100 


UO  14  1 - 1,10 
UO  14  J * 1,16 
15CLMT(I,J,rxABLB) 
CO & TIN U A 


SCALOB*SCxNl  41,0, 1‘IaoLE)  ♦ 0.  5 


11  (SCALOE  .Eg.  ICO.)  MBIT' E (0 , 16)  B BOIN.aSIN  (ITABLE) 
16  tOaii  A1  ( Iti-  , 10X  , • aUC  = ' ,F  6.2,6 A,  • SCALED  T.N111AL  SEBVaE 
1'  COAL  = ',  t5.2,12T,  'SCALaB  HE  AN  NA1X  IN  hUNLBEDTHS  •) 


IF  (aC*LGfl  .60.  1000.)  NEIXE  (6 ,20)  BEGIN  ,46 IN  (IX  ABLE) 
20  tCB6AX(lB-,  10  A,*  BEG  * ' ,F5.  2,bl,  ' SCALED  INITIAL  SEE  VAN 
1*  4.0AL  = ',F5. 2,1lX,':>CAi.ED  BEAN  BAIT  IN  1BGU5ANDTMS') 


MEIlt  (6, 

25  F0B6AT(1h 


551,'SCAt.ED  EtOCB  ' ) 


NE  II  a (6,30)  X VAL 

EG  Eh  AT  ( Iti  ,21X,S*5.  2,7f5.  1) 


wBlT  a (6,35) 

35  10 BEAT  ( 1H  ,211,16(11,' 


M All  a (6,40)  (ISCLNl  (1  ,J  ,ITABLN) 
4C  ruhtikl  (1b  , 1(1,  • mIENEB  ' , 1 6 ( 


NEITE  (6,45)  (ISCLhI  (2  ,J  , II  ABi.fi)  , J*1,1©) 
EoatiAX  (Id  , 101,'  6 Li  yUEUE  ',16(11,14)) 


UO  oC  1*3,6 

NEITE  (6,50)  AVAL  (I*  A),  (I6CT.NX  (I,  J,  11'ABt.E)  ,J*  1 , 16) 
foEfcAT(1h  , 101, *6  El  K*',t5.  a,1o(1I,14)) 

CONTINUE 


m Bat E (6,70)  (ISCL  mT  (10  ,J  ,IXAB  tJS)  , J*1,1o) 
tOftHAT  (IB  , 10X  , ' o ANNA  INPUT' ,1b  (IT, 14)  ) 


60  CONTINUE 


60  CONTINUE 


eltuen 


r r r r 


i' 


LiitlkU'ith  kuCiitbAci  FOB  *fkkbb Li  l’A rii-r.il,  £U0  i^UA*.  20  Out 


IdPiloil  hi-  AL*6  (A-n,  A,  L,  0-2)  , ikl  HiLk  (I,o,  a,  A) 

Luhttok  /lkuU'l/  1 if  AL  (loj/KVAi.  i7)  , iiAlX  (1C,1o,3)  , 21k  ( 3) 
CCtmOtt  /i. Ai?L AC/  V(Sb),(l,*l 
COBHC*  /ME*  J»AP/  10LkMt#dA2ri’k 

c 

h=b 

M*34 

KA2xlM=20 
lOLhbt*  1.i>-  2b 
C 

wO  bb  I PAGE =1,3 

hc.hu  it,  *0)  (2IM(xlABL£)  , IIAbb£=  1 , 3) 
k£Ab(5,20)  ('XVAL(J),  o=1,o) 
fcJLAb  p,2Q)  (‘XVAo(J),  J*7, 12) 

*EAi>45,zO)  (2 VAL  4 J)  , J»13,1b) 

20  tokflAi  (bi  10.C) 

C 

uO  50  iXAbLE* 1 ,3 
CALL  M2  a«k  (1'xAEL  b> 

CALL  AZBLl  (IX ABLE) 

CALL.  fixaEK  4llABi.ii) 

CALL.  LUG  AC.  (IXAoLk) 

50  COMlxMLE 
cA  Li.  kLQU'X 
bO  CCul'lMUE 
S'ioP 
EMD 


bitCCK  LAI  A 

LHLLLCll  bi.AL*6(A-b,K,L,0-Z)  , IMlkGEk  41  ,0,  cl,  k) 

CO KflO A /IMOUT/  TV  AL  (16)  ,K VAi.  47)  ,UAIT  (1C  ,1o,  3),Zlk  43) 
LA  I A K VAL/4  • ,3.  , x.  , 1. 5,  1.  , • S,  • 2/ 

£,M0 


LbbkoUlIkE  Lx*  V (P  ,M,I,FA  , V , A) 
ittfcLICll  k£  Ax.*fa  4 A- tt,U- 2)  , Ik'i&GEk  (i-M) 
CCabuM/LA  Pi.AC/£POCh  (10)  , APPaJA  (1C  ,0) 
uIMEUSlUM  t (25)  , F (25)  , U(2b),  V (oO) 

IF  (d.kwk)  GO  lu  50 
AH  * M/2 
£(1)  * 2 
L(2)  = 12 
F (1)  * 1 
*V)  = 1 
DO  1b  A « 2,  MH 

E (!♦  1).  * ( 2.  DC*  (x*I ♦1)*£(i)  ) / I 
F (!♦  1 ) » F (I)  *1 
10  COMT1MUE 
PUB  = Mb 
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U(1)  * 2/i  (MU) 

DC  2C  J * 2,  MU 

H (J)  * t(J»*P«S)  *£(*>))/£  (MH-J*1) 

*0  CO  MX  I MU  £ 

IS  lb*  * 2*  (Hu-  2*  ( MU/ 2)  > -1 
DO  40  I = 1,5 
V(I)  * 0 

JFIBSX  * (1*1)  /2 
JLAS'l  - alMO  (I ,Mu) 

DO  JO  J * JFI2SX,  DDAST 

V(l)  = vu)  ♦ U(uJ/4FU-Jt1>*F  (2*J-I*1)) 

3C  CO  MX  I MU  E 

V(I)  = 1j>IGM*V  (I) 

ISlGit  = ~i.SIGM 
40  cCmTIMUE 
fc  = M 

t>0  f A = 0 

A = . o :» 31 47  1 40559  945  30  54 17 232 121 4s UDO  / S 

LG  60  I*  1,  M 

FA  = f A ♦ V (I)  *P  (i*A) 

00  CGMX1MUE 
FA  = A*  FA 
mEXUuM 
£MD 


oU BkCUilUt  (llAoLE) 

IUPLICIl  «EAD*6(A-H,K,D,Q-Z)  , IM 1EGEB  (2,0,  tt,  M) 

CuHUOM  /IMG  UX/  X V AL  (1 6)  ,K VAi.  (7)  , WAIX  (1C, 16,  3)  , ZIM  (3) 
COttMuM  /i.  AFLAC/  V (50)  ,U,M 
CGEUCM  /ZO/  Z , OS'I  AuT 
EXlEfcUAL  PUNA 


2 - AIM  (IXABLE) 

DO  20  0*1,16 
X * 1VAD(J) 

CALL  LIM  V (P  MMU  ,N  ,X  ,XE£H3,  7,(1) 
UAIX  (1,  J,12ABLE)  * Z ♦ XEaU3 
20  Cu«XiMU£ 

BETUUM 

£ML 


DO  USD  £ P2ECISI0*  FUMC110M  PMMM  (S) 

IhtLICiX  *EAD*8  (A- fa,  K,  2,0-Z)  , IMXEGEB(I,J,tt,M) 
CCUHoM  /ZO/  Z,OS'IA£X 
«UEGA  * DSwfiX  (2*S) 

Phan  * DEAF  (-OHEG  A*Z)/ (d*OUEGA) 

«£XbaM 

£MD 


J 


r r. 


LU  o k CUT  I kE  PZMLl  (jAAcLE) 

IM  Pi*  xl  il  A£iAi«^d(&*hfhfi*i  0*i)  * Ik  1 bii  (I  pO  i A , ti) 

CUDHC«  /IkOUl/  X V Ai.  ( 1 o)  , tiV  A*.  (7)  , a AI'A  (10,1b,J)  , Elk  ( J) 

CGancfc  /u  APL  Ac/  V(bQ)#M,k 

CuMMCk  AO/  Z,OSIAfaI 

EXTEBfcAL  Pad,FMO  1,ufMDl 

OSIAEX  = 1 

* = dlk  (IlAu.Lt) 

tLZ  = uEXP(“l) 

iiO  lO  0*1,16 

X = TVAL(J) 

IF  (X  . LE  . a;  EZil  = 0. 

IF  (X  . GX.  Z)  CAi.1  L.tkV(PMd,k,'I,ttT,V,M) 

IF  (X  . Gl.  Z)  EZ1X  - UX  ♦ (I-Z)*LLZ 

BAIT  (Z,  J,I1  AaLE)  * Z ♦ EZIT 
20  coaxikOt 
*EXU«k 
END 


DOUElE  P&FLI£xOn  F UNCI  10  k Pttd  (a) 

itifLi C I X iih  Al^b  (k* H f K |LfO* / ik2'£^i£b  J|  b|  k) 

CO  ft  Ml  N /z&/  Z,OSXAtT 

coaaok  oj 

EX'iEZkA L FOHLl  ,DFML 1 
SS=  a 

CA  i.L  ZE2Q(F  CMLl,  DFMDl  , OuTAiii.  ,OMEGA) 

CSlAttT  * OMEGA 

Paul  = LEAP  <-GMEGA*Z)/ (E*QMEGA)  - DEX?  (-<.»(  1 *S)  )/ (S*a) 

fcEIUak 

END 

DOUc.lZ  pEEcIEIOk  FUNCTION  PMUl  (E) 

IMPLICIT  h£,Ai.*d(ii-fa,K,L,0-Zj  , IkiZGifi(I,.i,ii,k) 

ccnack  /zo/  z,osxaet 
Fad  = i>  - 1.  ♦ uexp(-S) 

BE'iUdk 

END 

DOUdLE  PdECIblON  F UNCTION  UFaUl(E) 

IMPLICIT  B£AL*8  (A-b,  A,L,Q-Z)  , INTEGEE  (I*u,  A ,-k) 

coaaok  /zo/  *,oexaat 

DFMDl  = l.-DEXP(-S) 

EETUaN 
Ek  L 

DCULtE  PRECISION  FUNCTION  FuMDl  (OMEGA) 

IMPLICIT  kEAL*8(A-h,K,  L,0-Z)  , INI EGEk  (I,G,M,k) 

COMMOk  SE 

FOMdI  * FMd  (OMEGA)  - EE 

aEXUkk 

Ek  D 
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sUBaOUliAt  lUaEA(IlAELZ) 

IBELICj.%  BEAL*8lA-H,K ,1,0-Z),  LulEOEE  (1,0, a,*) 

C Uhau  /lit  OUT/  1VAL<16),KVAL(7)  , WAIT  (10,16,3)  ,ZI»  (3) 
CuaBoa  /LAP LAC/  V(E0),a,N 
COMBck  /ZOK/  Z,OalA£T,K 
aXlEEEAL  £>6Eii,*BEK,D*HEK 

00  JC  1*3,  it 
h » AVAL  (1-2) 

OSTART  * 1 
A * ZxA  (IXkbLh) 

ELZ  * OEXPl-Z*(K+1)/&) 

00  20  J»1# 1 6 
’i  * T VAL  (J) 

If  (1  .LA.  X)  EZI1  * 0. 

If  (1  .01.  Z)  CAi.L  4.IMV(PatK,A  „C,Hl,V,a) 

If  (1  .01.  Z)  EZ1T  * ttl  ♦ (T-Z)*EL Z 
n All  (I . u ,11  ABLE)  * Z ♦ £ZI1 
20  CObTiAUE 
30  COaXIMUZ 
uETUa* 
chD 


DOUBLE  PaECIElO*  FUMCTIOA  PEEK  (S) 

IMPLItil  iBAL*8(  A-B/K/L/O-Z)  , 1M1EGEE  (1,0,  B,  h) 

COBBCE  /ZOK/  Z , OB1AET.K 
CCbBOa  BE 

EXIEaE AL  FOBEK ,0k  BEE 
EE  * E 

cAIL  ZtaC  (F  CBEK,  DFBEK, OblA*!  ,OBEGk) 

0SAAA1  - OBEOA 

PBEa  * LEAP  (-0HZGA*Z)/ (S*OB£OA)  - DEXP(-Z*  (S*  (A»1)  /A)  ) / (S*S) 

EElUfib 

E*0 


OuUBlE  PAECIEIOB  FUBCTIOB  FBEA(S) 

iaPi.!CIT  aEAL*8(A-d,K,*,0-Z)  , lal  EOEa  (1,0, K, a) 

CO  ABO*  /ZOK/  Z,0E1A&T,K 

FB EK  * E- {At1>*  (1- UA*1)/(A*nSi  ) **A) /K 

EElUAtt 

END 


lOUBlB  PaEClSIOM  FOACXIOli  OF  BE  K (E) 

IBPLlCiX  aEAL*8(A-ti,a,L,0-z).,  1BTEGEB  (1,J,B,B) 

COfldob  /ZOK/  Z ,OSXkkl , K 

OFdEA  * 1-((At1)/  (A*1zE))**(A»1) 

EETUdM 

EDO 


113 


n c 


cozaoA  si 

tCttEA  - fatMUBEoA)  - Si 
fiEID  aA 
LAD 
c 
c 

SUbaoUXIAE  AZGA8  (ilAbLE) 

latij-Cll  Ai,A*.*0(A-h,K#L#O-^  # laiEGsa  <i,z,  a,  *; 

COttflOA  /iAQUl/  l'V  Ai.  (Ibi  ,KVAi.  (7)  , AA1I  (10  ,1b  ,S)  ,ZIA  (3) 

CGE8GA  /LAP ZAC/  V4i0),8,il 

COCaCA  /<.0/  Z,OiTABl 

tXlEaAAZ  PGAa,FGAa,i>FGAA 

CilAEl  * 1 

i * AIV(IIASLE) 

C 

bG  20  J=1, 1 6 
1 = XVAi(J) 

If  (X  «EE«  EmxI  = C. 

It  (X  . UX.  Z)  CALi.  LIE  V (PGA  A #A  ,i  #EZITf  V 
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APPENDIX  C 

TABLES  OF  SCALED  EXPECTED  SERVER  LOAD 

All  tables  for  a specific  value  of  p are  grouped  together. 
The  tables  are  arranged  in  the  order  shown  in  this  list. 

Each  page  contains  results  for  three  values  of  z*. 


P 


2.0 


1.5 


1.4 


1.3 


1.2 


There  are  two  pages 
of  tables  for  each 
of  these  five  values 
of  p. 


z* 


1.0,  0.8,  0.6 

0.4,  0.2,  0.0 


1.1 

1.0 

0.9  There  are  three 

q g pages  of  tables 

for  each  of  these 
seven  values  of  p. 

0.6 

0.5 


4.0,  3.0,  2.0 

1.0,  0.8,  0.6 

0.4,  0.2,  0.0 
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2.00  SCALED  INITIAL  SERVER  LOAD  » 1*00  SCALEO  MEAN  WAIT  IN  THOUSANDTHS 

SCALED  EPOCH 

0*6)  0*80  1.00  1.20  1.40  1.60  1.80  2.00  2.50  3.0  3.5  4.0  S.O  6.0  T.O  8.0 
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1.50  SCALED  INITIAL  SERVER  LOAD  = 1.00  SCALEO  MEAN  WAIT  IN  THOUSANDTHS 

SCALED  EPOCH 

0.30  0.40  0.50  0.60  0.70  0.80  1.00  1.20  1.40  1.6  1.6  2.0  2.5  3.0  S. 0 8.0 
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This  research  provides  numerical  results  for  time-dependent  expected 
server  load  (mean  virtual  waiting  time)  in  single-server  queues  with 
Poisson  arrivals  and  gamma  distributed  service  times.  The  results  are 
presented  In  tabular  form  to  facilitate  their  use  by  practitioners  Involved 
In  the  study  of  operating  systems, 
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ABSTRACT  (continued) 

The  server  load  process  is  expressed  as  a net  input  process  (having 
stationary,  independent  increments)  modified  by  a reflecting  barrier  at 
the  origin.  In  queueing  systems  the  net  input  process  is  compound  Poisson. 
Two  other  choices  of  net  input  processes,  the  Wiener  process  (or  Brownian 
mStfoni  and  The  gamma  process,  provide  approximations  for  the  queueing 
process.  :>This  research  considers  groups  of  server  load  processes  whose 
parameters  are  selected  so  that  the  first  and  second  moments  of  their  net 
input  processes  are  matched.  An  existing  Laplace  transform  expression  is 
employed  to  obtain  transient  expected  server  load  at  specified  epochs.^ 

The  Laplace  transform  is  inverted  numerically  using  an  existing 
technique.  The  inversion  technique  is  first  applied  to  test  functions 
whose  exact  inverses  are  known,  and  the  results  are  also  compared  with 
queueing  results  obtained  by  other  methods.  These  investigations  indicate 
that  the  numerical  results  for  expected  server  load  have  four  to  six 
significant  digits  when  the  approximation  is  based  upon  thirty  four  values 
of  the  Laplace  transform  function.  The  computer  programs  are  coded  in 
FORTRAN  IV  using  extended  double  precision,  and  complete  documentation  is 
included. 

The  numerical  results  for  expected  server  load  are  tabulated  in 
scaled  form.  Each  of  the  93  tables  includes  results  for  a specific  combina- 
tion of  traffic  intensity  parameter  (twelve  values,  ranging  from  .5  through 
1.0  up  to  2.0)  and  initial  server  load  (either  six  or  nine  values,  ranging 
from  0 to  4 in  scaled  form).  An  individual  tables  has  results  for  the 
Wiener  process,  the  gamma  input  process,  and  eight  queues  with  different 
service  time  distributions;  each  of  the  ten  processes  is  evaluated  at 
sixteen  epochs.  Step-by-step  procedures  for  using  the  tables  are  included, 
and  several  sample  problems  are  presented. 

The  tabulated  results  allow  a comprehensive  study  of  the  error 
associated  with  using  the  Wiener  process  as  an  approximation  of  server  load 
In  queues.  This  study  confirms  that  the  Wiener  process  is  always  an  upper 
bound  and  that  the  approximation  is  best  for  queues  with  a traffic  intensity 
parameter  near  unity.  The  scaled  results  also  indicate  that  the  gamma 
input  process  and  queueing  process  with  deterministic  service  times  provide 
tight  lower  and  upper  bounds,  respectively,  for  expected  server  load  in  all 
queues  with  Poisson  arrivals  and  gamma  distributed  service  times. 
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